Correlation of CMB with large-scale structure: II. Weak lensing 
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We investigate the correlation of gravitational lensing of the cosmic microwave background (CMB) 
with several tracers of large-scale structure, including luminous red galaxies (LRGs), quasars, and 
radio sources. The lensing field is reconstructed based on the CMB maps from the Wilkinson 
Microwave Anisotropy Probe (WMAP) satellite; the LRGs and quasars are observed by the Sloan 
Digital Sky Survey (SDSS); and the radio sources are observed in the NRAO VLA Sky Survey 
(NVSS). Combining all three large-scale structure samples, we find evidence for a positive cross- 
correlation at the 2.5(T level (1.8a for the SDSS samples and 2.1a" for NVSS); the cross-correlation 
amplitude is 1.06 ±0.42 times that expected for the WMAP cosmological parameters. Our analysis 
extends other recent analyses in that we carefully determine bias weighted redshift distribution of 
the sources, which is needed for a meaningful cosmological interpretation of the detected signal. 
We investigate contamination of the signal by Galactic emission, extragalactic radio and infrared 
sources, thermal and kinetic Sunyaev-Zel'dovich effects, and the Rees-Sciama effect, and find all of 
them to be negligible. 

PACS numbers: 98.80.Es, 98.62.Sb, 98.70.Vc 



I. INTRODUCTION 



Great progress has been made in cosmology in the 
past several years, in large part due to measurements 
of the anisotropies of the cosmic microwave back- 
ground (CMB). Most recently, the Wilkinson Microwave 
Anisotropy Probe (WMAP) experiment has generated 
high-resolution, all-sky maps of the CMB P, Q . While 
a great deal of attention has been given to the "pri- 
mary" anisotropies in the CMB imprinted at z ~ 10"^, 
experiments such as WMAP are also sensitive to sec- 
ondary anisotropies caused by electron scattering or grav- 
itational potentials at lower redshifts. Among these 
secondary effects are the thermal and kinetic Sunyaev- 
Zel'dovich (tSZ/kSZ) effects, the integrated Sachs- Wolfe 
(ISW) effect, and gravitational lensing. These secondary 
anisotropies are important for two reasons: first, they can 
act as "foregrounds" for primary CMB studies if they 
are not adequately modeled; and second, they provide 
information about the growth of structure at low red- 
shift. Since the secondary anisotropies are subdominant 
on the large angular scales observed by WMAP, they are 
most easily detected by cross-correlation with large scale 
structure (LSS). Several groups have cross-correlated 
WMAP data with LSS tracers in order to study the 
ISW and tSZ secondaries and extragalactic point sources 
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Lensing of the CMB is a secondary anisotropy that 
has attracted considerable attention because of its po- 
tential to probe the matter distribution at intermediate 
redshifts z ~ 0(1) [2l|, [H, [H . Lensing of the CMB 
has several potential advantages as a cosmological probe 
relative to using lensing of galaxies: (1) it probes a red- 
shift range inaccessible to galaxy lensing surveys since its 
"source screen" is at z sa 1100; (2) its redshift is known 
to high accuracy; and (3) as a Gaussian random field the 
CMB does not suffer from intrinsic alignments. Lens- 
ing of the CMB does however have its limitations: (1) 
the low signal-to-noise ratio with current data; (2) one 
cannot make tomographic measurements of the redshift 
dependence of the lensing signal without including exter- 
nal data sets [23|; (3) foregrounds; and (4) instrumen- 
tal systematics are a problem, especially for the lensing 
auto-power spectrum. 

In the present paper, we aim to measure the weak grav- 
itational lensing effect in the WMAP data. We do this 
in cross-correlation with large-scale structure (LSS) in 
order to improve the signal-to- noise ratio. We recon- 
struct the lensing deflection fie ld |ll2| | from WMAP using 
quadratic estimator methods '25*1. We then measure the 
cross-power spectrum of this deflection field with lumi- 
nous red galaxy (LRG) and quasar maps obtained from 
the Sloan Digital Sky Survey (SDSS) and with the radio 
source maps from the NRAO VLA Sky Survey (NVSS). 
The detailed sample definitions, an analysis of the ISW 
effect, and parameter constraints are presented in a com- 
panion paper ( "Paper I" [2^ ) . 

There have been two previous searches for the len sing 
effect on the CMB in cross-correlation. Hirata et al. [23 
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used the the first-year WMAP data release and a smaller 
sample of LRGs from the SDSS. They found a nondetec- 
tion: 1.0 ± 1.1 times the expected signal. More recently, 
Smith et al. [Ill used the third-year WMAP data and 
the NVSS radio sources (although with different selec- 
tion cuts than those used here); they find a 3.4(7 result 
(1.15 ±0.34 times the expected signal, although their ex- 
pected signal is 15% larger due to a different bias, redshift 
distribution, and fiducial cosmology). The analysis pre- 
sented here draws heavily on that of Hirata et al. [27| . 
but benefits from: (i) reduced instrument noise due to 
the longer integration time in the 3-year WMAP data 
release; (ii) a larger sky area available in the SDSS; (iii) 
the inclusion of SDSS quasars and NVSS sources in ad- 
dition to the LRGs; and (iv) improved treatment of tSZ 
and point source contamination. The current analysis 
also extends these previous analyses in that we put a lot 
of effort in accurate determination of the bias weighted 
redshift distribution, using cross-correlation information 
with other galaxy samples. Even though this is not re- 
quired for claiming a detection, it is needed for any cos- 
mological interpretation of the signal and thus for any 
consistency check of cosmology based on this analysis. 

The analysis of this paper will be mostly based on a 
"fiducial" cosmology, which is the 3-year WMAP-only 
best-fit 6-parameter cosmology 29] from the point in the 
WMAP Markov chain with the highest likelihood. The 
basic objective is to measure the strength of the lensing 
signal, determine its statistical error and sensitivity to 
systematics, and establish whether it is consistent with 
expectations for the fiducial cosmology. We explore the 
effects of cosmological parameters on the lensing signal 
by running a Markov chain in Paper I. The only place 
where we vary cosmological parameters is in some of the 
extragalactic foreground tests to show that even extreme 
variations in the assumed cosmology do not affect the 
basic conclusion that the foregrounds are negligible. The 
parameter sets we will use are as follows. The "fidu- 
cial" (WMAP) cosmology is fiat with fJ^/i^ = 0.0222, 
Vtmh'^ = 0.1275, h = 0.727, CTg = 0.743, and = 0.948. 
The "high-f2m" is an extreme model that maintains the 
angular diameter distance to the CMB and keeps O&Zi^, 
flmh^ , n-s, and the primordial power spectrum the same 
but increases f2,„ until there is no cosmological constant 
(f^A = 0, n„i = 1.33, h = 0.31). The "high-erg" model 
is the WMAP + SDSS vaniUa model from Tegmark et al. 
^ and has Qbh^ = 0.0232, n„,h^ = 0.1454, h = 0.695, 
as ~ 0.917, and Ug — 0.977. The last model would 
give a large tSZ effect since it depends strongly on the 
normalization. Of these models we consider the fidu- 
cial and high-(T8 models to be viable, while the high-fim 
model is already strongly ruled out by supernovae, large- 
scale structure, and measurements of the Hubble con- 
stant [S^, HH [13, HI] and is included only to make the 
point that the foreground analysis is insensitive to cos- 
mology. 

The outline of this paper is as follows: the theoretical 
predictions for lensing of the CMB are reviewed in Sec- 



tion HH The WMAP and SDSS data, including the LSS 
samples used, are described in Section IIIII The analy- 
sis methodology is presented in Section |TVl We present 
our results in Section |V] and systematic error estimates 
in Section Extragalactic foregrounds deserve special 
consideration and occupy Section IVIII We conclude in 
Section IVIIII The appendices cover the description of 
point source contamination in bispectrum language (Ap- 
pendix |X| ; the halo model description of the bispectrum 
(Appendix |B]); the kSZ and Rees-Sciama foregrounds 
(Appendix \Q ; cross-correlations of different foreground 
components (Appendix |D]); and the weak lensing likeli- 
hood function (Appendix [E]) . 

II. THEORY 

Gravitational lensing re-maps the primary CMB signal 
according to the relation 

T(n) =f(n + d(n)) (1) 

Here T(n) is the observed CMB temperature in direc- 
tion fi, T is the unlensed (primary) CMB temperature, 
and d is the lensing deflection angle. In the case where 
the lensing is adequately described by a single deflection 
(the "Born approximation"), one may deflne a lensing 
potential <& and convergence k satisfying d — V<i> and 
K — — V^$/2, either one of which contains the full infor- 
mation in the lensing fleld [3^ ] . The convergence is given 
by 

nXCMB 

k = 4:TtGpq / (H-z)sinl-x(cotKX-cotK XCMB)rfx> 
Jo 

(2) 

where po is the mean density of the universe today, x is 
the comoving radial distance, and sihk X ^^'^ coIk X are 
the sinelike and cotangentlike functions (x and 1/x in a 
flat universe). The correlation of the galaxy density with 
the convergence field can be determined by the Limber 
equation, which yields 

AnGpo f dN 
^ J{dN/dx)dx J ^ - 

x{l + z)Pik = l/x)dx (3) 

[27| . where dN/dx is the comoving distance distribution 
of the galaxies, P{k) is the matter power spectrum, and 
bg is the galaxy bias. To measure the cross-spectrum , 
we need to be able to reconstruct k from the WMAP 
data. If we knew the primary CMB signal T(n), this 
would be a simple exercise. In the absence of this knowl- 
edge, we must rely on its statistical properties in order to 
do a lensing analysis. The primary CMB is a Gaussian 
random field with power spectrum Ci, i.e. 

{'^hmi'^hm2) = C^h^hh^mim^: (4) 

where Ti-^mi is a multipole moment of the primary tem- 
perature field. Expanding Eq. ([1]) to first order in d. 
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and working in harmonic space, one can show that the 
corresponding two-mode expectation value for the lensed 
CMB is 

LM 

X ( I ^^LM■: (5) 

\ —mi 1712 M I ^ ' 

where the lensing coupling coefficient [131 is 



^ _ v/(2£ + l)(2^i + f)(2;2 + l) (h h L\ 
^"-'''^ LiL + l)V^ [O Oj 

x{[L{L + 1) + hih + 1) - hih + l)]Ci, 
+ [L{L + 1) - hih + 1) + hih + 1)]^}. (6) 

In Eq. (O, the first term is simply the unlensed expec- 
tation value, and the second term represents off-diagonal 
correlations (i.e. correlations between modes with dif- 
ferent I or m) induced by lensing. The lensing recon- 
struction techniques that we will describe in Section FlVAl 
are based on an optimal weighting of these off-diagonal 
terms. 



highly elliptical beams do not preserve high-Z informa- 
tion. The inclusion of these bands would only slightly re- 
duce the total (CMB-I- noise) power spectrum Ci+Ni: for 
3 bands (Q-f V-l-W) the isotropically averaged C/ -I- Ni is 
reduced by a maximum factor of 1.10 versus V-l-W only. 
Including the K and Ka bands as well would improve 
this factor to 1.13 (relative to V-fW). This maximum is 
reached at I ~ 400 since at lower I the sampling variance 
dominates and at higher I the K, Ka, and Q band beam 
transfer functions are too small. Since the lensing recon- 
struction is quadratic in temperature, the factor of 1.10 
in temperature variance translates into at best a factor 
of 1.10 reduction in the standard deviation of Cf (the 
improvement is even less than that if not all lensing in- 
formation comes from I ~ 400 where the Q band gives 
the most improvement). There is also the issue of point 
sources, which were already a major concern in our first- 
year analysis [27j and would presumably become worse 
as the 3-year WMAP data reduces the noise at high I: 
because the lensing reconstruction is quadratic in tem- 
perature, flat-spectrum radio sources would produce 4.4 
times as much contamination in Q band as in V and 19 
times as much in Q as in W. We therefore decided to use 
only V and W bands for the main lensing reconstruction, 
and use K, Ka, and Q only for foreground tests. 



III. DATA 

A. Cosmic microwave background from WMAP 

We use the first three years of CMB temperature data 
[U, 0, [M m from the WMAP satellite [11, [13 located 
at the Sun-Earth L2 Lagrange point. WMAP carries a 
set of ten differencing assemblies (DAs) that measure the 
difference in microwave intensity between two points on 
the sky. The satellite rotates through an interlocking 
scan pattern that allows each DA to build up a map of 
the entire microwave sky. The DAs are designated Kl, 
Kal, Ql, Q2, VI, V2, Wl, W2, W3, and W4, where the 
letters indicate the central frequency (K, Ka, Q, V, and 
W correspond to 23, 33, 41, 61, and 94 GHz respectively). 
A new sky map is produced by each DA every year; since 
this analysis uses three years of WMAP data, this means 
there are 30 sky maps available. 

The WMA P maps have been generated in the 
HEALPix |ll3l | pixelization system ^38| at resolution 9, 
which has 3 145 728 pixels each covering a solid angle 
of 47.2 arcmin^. The actual resolution achieved in the 
maps is determined by the beam and hence varies with 
frequency, with the lower frequencies giving poorer reso- 
lution. The lensing analysis relies on the high multipoles 
il > 400) that are only accessible to the V and W band 
DAs (which have 21 and 14 arcmin beams respectively), 
and hence we will almost exclusively use the maps from 
these bands. The K, Ka, and Q band maps are never 
fed through the lensing pipeline because their large and 



B. Large scale structure from SDSS 

Two of the large-scale structure samples used in 
this analysis are photometric luminous red galaxies and 
quasars selected from SDSS imaging. The SDSS drift- 
scans the sky in five bands iugriz) [3^ under photomet- 
ric conditions [43,[4l| using a 2.5-meter optical telescope 
[4^ with 3 degree field of view camera [i^ located in 
New Mexico, USA [4l|. The photometric and astromet- 
ric calibration of the SDSS and the quality assessment 
pipeline are described by Refs. [3, [ii, [12 [13, [11 , re- 
spectively. Bright galaxies [49l. luminous red galaxies 
[50|, and quasar candidates [5l(] are selected from the 
SDSS imaging data for spectroscopic follow-up with a 
spectrograph connected to the same telescope 52]. We 
only use the imaging data since the spectroscopic galaxy 
sample covers only the very low redshifts (which are inef- 
ficient for CMB lensing) and the number density for the 
spectroscopic quasars is too low. Additionally in angu- 
lar cross-correlation there is no advantage to having pre- 
cise redshifts. The selection criteria for our LSS samples 
are related to, but distinct from, those for spectroscopic 
target selection. The SDSS has had seven major data 
releases [H, [H, [H [M [H [M [H • 

The use of SDSS LRGs and quasars as LSS tracers 
was driven by a compromise among several competing 
requirements. First, there is a need for large sky cov- 
erage to improve statistics. Second, one desires a sam- 
ple with high number density and bias to minimize the 
degradation of the correlation coefficient due to Poisson 
noise. Third, one would like the LSS sample to cover 



4 



a broad range of redshifts, since the CMB is lensed by 
structures at all redshifts and use of an LSS sample with 
small Az implies a small correlation coefficient with the 
convergence field at z — 1100. 

The photometric LRGs observed by SDSS satisfy all 
three of our desirata; we have used the same selection 
criteria as were used in the previous lensing '27| and an- 
gular clustering [^] analyses, although the sky coverage 
has increased by a factor of 2. The photometric LRGs 
cover the redshift range 0.2 < z < 0.7, which is signif- 
icantly deeper than the SDSS spectroscopic galaxy and 
LRG samples. However most of the CMB lensing arises 
from structures at higher redshift. Therefore in this anal- 
ysis we have also used the SDSS photometric quasars. 
These can be seen from much larger distances than LRGs 
and cover the redshift range z < 2.7, which broadens the 
overall redshift coverage of our LSS samples. However 
their shot noise is greater: the LRGs have a comoving 
number density of 4 x 10~^/i^ Mpc~^, versus a maximum 
of 10~^h^ Mpc~^ for the quasars. Of course, maximum 
theoretical signal-to-noise (though not necessarily the ob- 
served signal-to-noise) is achieved by combining the LRG 
and quasar constraints, and including the NVSS sources 
as described in the next section. 

For this analysis we included the photometric SDSS 
data obtained between 1998 September 19 and 2005 June 
11, and gridded in HEALPix resolution 10 due to the ex- 
tensive small-scale structure in the survey boundaries. 
The survey area after rejecting regions of high reddening 
or stellar density, bad seeing, and regions contaminated 
by bright stars is 2 025 731 pixels (6641 deg^) for the 
LRGs and 1842 044 pixels (6039 deg^) for the quasars. 
The difference is due to the more stringent reddening cut 
E{B — V) < 0.05 used for the quasars, as compared to 
E{B - V) < 0.08 for the LRGs. 

The LRG and quasar samples in Paper I were sliced 
into thin redshift slices in order to measure the redshift 
evolution of the ISW effect. The CMB lensing window 
function is broad and so we do not require high redshift 
resolution, thus for this paper both LRG samples in Pa- 
per I were grouped into a single LRG catalog, and both 
quasar samples were grouped into a single quasar cata- 
log. The number densities are 129 LRGs and 39 quasars 
per square degree. The details of the selection criteria 
and determination of sky cuts can be found in Paper I 
and will not be repeated here. 



C. NVSS sources 

The NVSS radio sources are useful for the lensing anal- 
ysis because the survey covers a larger portion of the sky 
than SDSS (27 361 deg^) and is deeper than the LRGs: 
the median redshift is Zmcd ~ 1-0. The sample is iden- 
tical to that considered in Paper I and its detailed con- 
struction will not be repeated here. The main points are 
that the sources are selected to have flux > 2.5 mJy in 
L-band (1.4 GHz) and are unresolved at the Very Large 



Array (VLA) in the "D" configuration used for NVSS. 
(The catalog goes down to 2.0 mJy, however it is 50% 
complete at 2.5 mJy and there is a danger of spurious 
power if one goes to very low completeness levels without 
a detailed investigation of the noise properties of the sur- 
vey.) The survey covers the whole sky except for regions 
with S < —40°. We also removed heavily contaminated 
regions in the Galactic Plane (|&| < 10°) and regions 
contaminated by sidelobes from bright sources. The sur- 
vey details can be found in the technical paper by Con- 
don et al. [6l|. These cuts accept 1 104 983 sources. As 



with the ISW analysis, the NVSS sources are treated as 
a single slice here because one does not have meaningful 
photometric redshifts. Indeed, even the approximate de- 
termination of the redshift distribution was a significant 
challenge (see Paper I). 

The bias-weighted redshift distribution f{z) = 
bdN/dz for the NVSS sources is obtained by fitting a 
F-distribution to the cross-correlation of NVSS with the 
SDSS samples and with sources from the 2-Micron All 
Sky Survey (2MASS) [63|. The procedure is described 
in great detail in Paper I. Note that we did not use the 
NVSS auto-power in our fits as it appears to contain 
spurious power whose angular spectrum is not precisely 
known. 



IV. ANALYSIS 

A. CMB lensing reconstruction 

The most straightforward method to correlate LSS 
with lensing of the CMB is to construct a map of the 
lensing convergence from the CMB data, and then com- 
pute the cross-power spectrum with the LSS data. In 
practice the convergence map produced by existing re- 
construction techniques is a highly nonlocal function of 
the CMB data, which represents a problem when cuts 
due to point sources or the Galactic plane are taken into 
account. Therefore, we will pursue a slightly different 
strategy, namely to construct a vector field v from the 
CMB data that is correlated with the deflection field, 
but which only depends locally on the CMB tempera- 
ture (i.e. there is very small dependence on points more 
than a few degrees away). We apply Galactic and point 
source masks to the field v instead of directly to the tem- 
perature map. We then find the cross-power spectrum of 
V with each LSS sample by standard quadratic estima- 
tion techniques. The same basic concept was used in 
^2^, but some modifications have been made here so we 
describe the current version of the algorithm in detail. 
The methodology is outlined in Section IIV A 11 and the 
choice of weight functions is determined in Section llV A 21 
The final construction of the lensing maps, including the 
masks, is described in Section llV A 31 All real-space steps 
in the lensing reconstruction use HEALPix resolution 9. 
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1. Methodology 

We perform a reconstruction of the CMB lensing field 
using a quadratic estimator. Quadratic estimators have 
been proposed by many authors for lensing reconstruc- 
tion from the CMB temperature [H, HI, H, [Hi , CMB 
polarization 2ll. 66^ 671 68| . and diffuse high-redshift 21 
cm radiation 69|, [70|, [7l| . For the reconstruction based 



on CMB temperature anisotropics such as those observed 
by WMAP, quadratic estimators provide similar signal- 
to-noise to more complicated likelihood-based methods 
[t^ I- although this will not necessarily be the case for 
future CMB exp eriments that are sensitive to _B-mode 
polarization [73|. The three key steps in constructing a 
quadratic estimator, as described in [25| are as follows. 

1. The temperature field is filtered to produce a 
weighted temperature map, 



Wlra = 



Tim (observed) 



(7) 



where Bi is the beam transfer function (including 
both the physical beam and the Hcalpix pixel win- 
dow function) and Cf* is a weighting function. In 
the case considered by [IE] where the noise is statis- 
tically isotropic, it is optimal to use C;^* — Ci+Ni 
where Ni is the noise power spectrum; as discussed 
below statistically isotropic noise is not an appro- 
priate assumption for WMAP. 



2. The filtered temperature gradient is produced: 
G(n) = V^C';VFi™yi™(n), 



(8) 



where Ci is the unlensed CMB power spectrum (or 
one's best estimate of it). 



3. A "temperature- weighted gradient" 
G(n) = M^(fi)G(n) 



(9) 



is computed, and in the methodology of 25| this is 
filtered to produce a deflection angle or convergence 
map. 

One can show that the end product of these manipula- 
tions is 



l'l"ra'm" 

I V I" \ Ti,m,{ohs)Ti„m"{ohs) 
BvBinCJi^Cp 



.(10) 



In the absence of lensing [n — 0), this will have ex- 
pectation value zero because for / ^ we either have 
V — I" and m' = —to" (in which case either / is even 
and the 3-j symbol vanishes, or / is odd and Jwi" = 0), 
or else (Ti/™' (obs)Ti//m// (obs)) = 0. This is essentially 



a consequence of statistical isotropy. A particular re- 
alization of the lensing field breaks statistical isotropy 
and induces off-diagonal terms in the covariance matrix 
Eq. ([5]). These add coherently to yield a contribution to 
{Girn), at linear order in k this will be proportional to 
Kim by rotational symmetry. 

In practice, this methodology runs into several prob- 
lems, and we make several changes in order to avoid 
them. One problem is the leakage of the bright fore- 
grounds in the Galactic plane into the survey region when 
Eqs. ([7]) and ^ are implemented. We solve this by 
setting the portion of the temperature map within the 
Galactic Kp2 cut to zero before computing Eq. This 
results in artifacts near the Kp2 boundary but as we will 
see most of the sky is uncontaminatcd. 

The second problem is that for non-uniform noise such 
as that in WMAP, the above procedure turns out to be bi- 
ased because the same instrument noise appears in W^(n) 
and in G(n), leading to a bias in G(n). [We can also see 
this from Eq. pop . since T^m' and Tii'm" contain noise 
correlations if the noise varies across the sky, which is 
certainly the case for WMAP.] If the noise covariance ma- 
trix is known well this bias can be subtracted off. The 
WMAP noise properties are clean enough that this is 
probably possible (especially if the resulting lensing map 
is to be used in cross-correlation). However the simplest 
and most robust way to avoid this problem without being 
sensitive to possible errors in the noise covariance matrix 
is to use cross-correlations among the 18 V and W band 
maps obtained by WMAP. Letting Greek indices a, (3, 
etc. denote the maps, we may replace Eq. ^ with 

G"'3(n) = i[M/"(n)G'^(n) + M^'5(fi)G"(fi)], (11) 

where W"'{n) represents the map produced by feed- 
ing the ath WMAP map into the pipeline that con- 
structs W{n). This procedure provides 18 x 17/2 = 153 
temperature-weighted gradient maps. Since we have 18 
maps there is little loss in using cross-correlations only 
(e.g. if the maps all had the same noise then the loss of 
signal-to-noise ratio ranges from 1 in the sampling vari- 
ance limited case to ^18/17 in the noise limited case); 
this is a small price to pay for eliminating the reliance on 
the detailed noise covariance. 

The third problem with the standard quadratic estima- 
tor that comes up is that the filtering of G(n) to obtain 
an unbiased estimator of k is highly nonlocal. This is an 
issue because it spreads artifacts from point sources over 
the entire sky. Our solution to this is to instead apply a 
Gaussian filter to G(n) to obtain a filtered vector field 
v(n): 



"im ~ ^ '^Im 



(12) 



Here || and _L represent the longitudinal (vector) and 
transverse (axial) multipoles, which are spin-1 analogues 
of the E and B multipoles for tensor fields. The filtering 
scale tJo is set to 10"^ radians (34 arcmin). The fields 
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Response of lensing estimator 
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FIG. 1: The response factor Ri of Eq. ([T5j. This is the cal- 
ibration relating the (mean of the) reconstructed map v to 
the underlying convergence factor k. Note that it depends on 
scale. 

(n) are suitable for cross-correlation studies since ar- 
tifacts from point sources are local and can be masked. 
Once they are constructed, the v"^(n) are averaged to- 
gether to make a final lensing map according to 

v(n) = . (13) 

(The weights Wap will be specified in Section HV A 21 ) 

The use of the field v instead of a simple unbiased 
estimator for k means that some care must be taken 
in intepreting the LSS-v cross-power spectrum. Specif- 
ically, V is an unbiased estimator for some filtered ver- 
sion of K rather than k itself, and this filtering must be 
accounted for to obtain meaningful results. Given a par- 
ticular convergence field, the multipole moments of v"'' 
can be shown to have expectation value 

KH^^^^) = Rl^lm, (14) 

where the response factor Ri is 



(cf. Eq. 17 of Ref. ^^). This can be derived by plug- 
ging the expectation value from Eq. ([5]) into Eq. PH)) . 
using the 3-j symbol orthogonality relations to collapse 
the sums over m' and m", and incorporating the Gaus- 
sian factor from Eq. This response factor is plotted 
in Fig. m 

2. Weighting 

Finally, we note that the implementation of Eqs. ([7]) 
and ([5]) requires an estimate of the unlensed CMB power 



spectrum C; and a choice of weighting function Cf^ . For 
Ci we use the CMBFast [tII prediction for the fidu- 
cial cosmology. The "actual" C; estimates from WMAP 
could have been used instead, but this would have the un- 
desirable feature of introducing a long-range dependence 
of G(n) on the CMB temperature, which is a problem 
if the Galactic plane has been removed. We note that 
in any case an error in the assumed Ci does not lead to 
any bias in the convergence field; rather it changes the 
response factor Ri This is equivalent to a (possi- 

bly scale-dependent) multiplicative or calibration bias in 
the lensing reconstruction, which cannot give a spurious 
signal. 

For Cj**, we use an approximation to the total power 
spectrum (CMB-|-noise) obtained as follows. Each of the 
18 maps has a noise variance per pixel tr^ j that depends 
on both the map and the pixel i. In the limit where 
the noise is uniform and we use the Hu version of 
the quadratic estimator, it would be optimal to define 
Cf* :^Ci + Ni. We have thus chosen to write 

Q"* = Ci + (0.015655 AiK2)(B,vi)-2(i _ f/uoo^)-^ 

(16) 

for I < 1200 and oo for I > 1200. The second term in the 
I < 1200 case is a good approximation to the effective 
noise in the I « 400 range where Ni « C/. Equation 
p6| has the advantage of having a cutoff in CJ'* at 
high I, so Eq. ([7]) can be computed by multiplication in 
harmonic space, and this cutoff is smooth, which avoids 
ringing when converted back to real space. 

Finally we come to the selection of the weights Waf3- 
If the WMAP beams were all the same and we were in- 
cluding the auto-correlations as well as cross-correlations 
in the lensing reconstruction - i.e. if we included a = (3 
terms in Eq. I|13p - then the optimal weight would be 
simply Wap = \/[NaNp), where Na represents the noise 
variance in map a. For simplicity we choose to use this 
weighting scheme, where for we have used the noise 
at the I ~ 400 multipole (i.e. we take the mean-square 
noise per pixel and multiplied by the inverse-square beam 

-^400 ) • 

3. Masks and final lensing maps 

The lensing maps v have two major types of arti- 
facts: one in the Galactic plane, induced by the Kp2 
cut applied in Section IIV A II and by foreground emis- 
sion; and spurious features surrounding the point sources. 
Both of these must be masked before proceeding. The 
Galactic emission is treated by masking out all pix- 
els within the KpO region, or within 5 degrees of its 
boundary. This "Kp05" cut was also used in our first- 
year analysis [23| and accepts 2 064181 HEALPix res- 
olution 9 pixels (/sky = 0.656). We also remove any 
pixels centered within 2 degrees of a WMAP-detected 
source, as determined in the third- year catalog [T!] , which 
gives a final "Kp052" mask that accepts 1 825 036 pixels 
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FIG. 2; The CMB lensing map V ■ v'^^' in Galactic 
MoUeweide projection, degraded to 40 arcmin resolution. The 
gray region is rejected by the Kp052 mask. The scale runs 
from -2 X 10* (black) to +2 x 10** (white). The map is domi- 
nated by noise (both instrumental and due to the finite num- 
ber of primary CMB modes available for the reconstruction) . 



(/sky = 0.580). Our main results will be shown with the 
Kp052 mask. 

In order to study the frequency dependence of the sig- 
nal, we construct individual- frequency maps v^^, v^'^, 
and v^^, which are obtained by applying Eq. but 
restricting the averaging to products of maps in the spec- 
ified bands. For our weighting, we find that VV, VW, 
and WW contribute 26.35%, 53.24%, and 20.41% respec- 
tively, i.e. 

v(n) = 0.2635v^^(n)-h0.5324v^'^(fi)-t-0.2041v'^^(n). 

(17) 

The reconstructed lensing map is shown in Fig. [51 



B. Cross-correlation 

1. Cross-spectrum estimation method 

We have cross-correlated the lensing maps v with 
the LRG and quasar maps X using a quadratic cross- 
spectrum estimator. These methods have been exten- 
sivelydeveloped for both LSS and CMB applications 
[Zi,lli,[z3|- We write the LSS map as a vector of length 
■^pix,LSSj and the CMB lensing map as a vector of length 
2-/Vpix,CMB (since it has two components). The cross- 
spectrum is written in the form of a matrix, 

c'if = {g^v,) = Y.Ri^rYCmmL{j). (is) 

Im 

Here C^^^ is an A'^pix^LSS x 27Vpix,CMB matrix; Hs a pixel 
index in the LSS map; and j encodes both a pixel index 
in the CMB lensing map and a unit vector on the celestial 
sphere {eg or e^). In the second equality we have written 
as shorthand Yim{i) for the value of a scalar spherical har- 
monic at pixel i, and i^l'„(j) for the component of a vec- 
tor spherical harmonic at the pixel and in the direction 
corresponding to index j. This matrix is assumed to be 



a linear combination of templates, C\j = J^A'^^^Aij, 
where the templates Paij are known and their amplitudes 
are to be estimated. We then construct the quadratic 
estimators 

= g^w(f)PAw(^V (19) 

and 

^^^s = Tr[w(s)PAw(^)pi;], (20) 

where w*^^-* and w'^'^^ are symmetric weight matrices 
(A^pixiss X A^pix,LSS and 2A^pix,CMB x 2A^pix,CMB respec- 
tively) . 

The theory of quadratic estimation [77| provides two 
key results concerning Eq. (|5D)) . The first is that the 
quantities c"^ = [F~^]'^^qB are unbiased estimators of 
the C^. Note that this is true regardless of the (non- 
)Gaussianity or correlations between g and v, and re- 
gardless of the choice of weight matrices (so long as F 
is invertible) . The second key result concerns optimality. 
If g and V are Gaussian and weakly correlated, then the 
quadratic estimator is minimum- variance with the choice 
of weight matrices w's) = [C'-s)]-^ and w(^) = [C(^)]-i. 
In our particular application the assumption of Gaus- 
sianity is violated, because v is constructed from prod- 
ucts of the CMB temperature field and hence is itself 
non-Gaussian. However, in the absence of better infor- 
mation, we use this optimality result to guide our choice 
of estimator. The actual variance of the estimator will 
be given not by the matrix F (which in our context need 
not be the Fisher matrix) but rather by simulations (see 
Section HVCl). 

We have used two sets of templates P in the cross- 
correlation analysis. The first set of templates uses band 
powers, i.e. 

pm, = E E RinnmLij)^ (21) 

in which case the coefficients can be interpreted as 
the values of Cf^ in the range Zmin(^) < I < ^max(^)- 
We use 13 bands covering the range from 2 < I < 400. 
The other method is to use as a template the theoretical 
signal for the fiducial cosmology. 

2. Weights 

We require a weight matrix for the LRGs, the quasars, 
the NVSS sources, and the CMB lensing map v. The pre- 
scriptions below for the LRGs and quasars are believed 
to be a good approximation to the inverse-covariance ma- 
trix. For the CMB lensing map, the true covariance ma- 
trix is not known and w*^'^) should be thought of only as 
a heuristic weighting scheme. 

For the LRGs, quasars, and NVSS sources, a prelimi- 
nary power spectrum was measured in Paper I. We have 
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thus estimated the LRG covariance matrix as 

cr^^ = ^+ E cr^'^Yr^muji (22) 

Im: i<400 

where the first term is the Poisson noise contribution, and 
second term is the angular clustering. The weight matrix 
is then w(^^^) = [C^^^'^)]"^ The Poisson noise depends 
on the mean number of LRGs per pixel, n = 0.424. Note 
that n is not allowed to depend on the pixel number, as 
this would place more weight on overdense regions and 
thus potentially bias the results. The clustering contri- 
bution to the weight matrix is cut off at ^ > 400, where 
Poisson noise dominates the LRG power spectrum, and 
we set Co — Ci = 10~^ in order to reject the monopole 
and dipole. A similar method was applied to the quasars 
and NVSS sources. 

For the CMB lensing map, we followed a similar pre- 
scription for constructing a covariance matrix C*-^^ using 
the power spectrum of v taken from a simulation. How- 
ever there are three differences that appear when han- 
dling V instead of an LSS map. First, v is a vector so 
there are both longitudinal and transverse power spec- 
tra. Secondly, the power spectrum of v is in fact variable 
over the sky, with higher values at low ecliptic latitude 
where the WMAP maps are more noisy. Thirdly, v does 
not have "shot noise" on small scales, which makes the 
covariance matrix C''*'-' nearly singular. We can work 
around all of these problems by writing 

1 771 

E[C;^^'^ - (23) 

Im 

The use of both longitudinal (||) and transverse (_L) 
modes is necessary since v is a vector. The factors ipi 
attempt to weight different parts of the sky depending 
on their noise. Finally, the first term has been intro- 
duced to ensure that C has no eigenvalues less than A/'/fi, 
where the noise floor is A/" = 2 x 10^ and the pixel area 
is 57 = 4.0 X 10^^. The sum over I runs only over those 
values for which cf"^^ > Af; this is / < 265 for the longi- 
tudinal modes and / < 267 for the transverse modes. (In 
these cases we take out the noise floor JV.) The CMB 

lensing weight is then w^"*'^ = [Cipprox]"^. 

We also need to select the factors tpi, which control the 
relative weighting of different regions of sky. Since ideally 
we would like Eq. (|23p to be the true covariance matrix 
of V, it is desirable to have i/jf proportional to the local 
power spectrum of v, which is highest near the Ecliptic. 
We therefore divide the sky into several regions based on 
ecliptic latitude /3, and computed pseudo-C; power spec- 
tra of V (again taken from a simulation rather than the 
real data in order to guard against possible correlations 
of the data with the weight matrix) in each. It is found 
that these power spectra are ~ 25% higher than the av- 
erage near the Ecliptic, and ^ 25% lower than average 



near the ecliptic poles. Therefore we choose the heuristic 
weighting factors 

= Vl-25- 0.50|sin/3,|. (24) 

We have not chosen to do a more sophisticated weight- 
ing (e.g. more general dependence on /3, or allowing 
the weighting factors to depend on /), noting that the 
quadratic estimator is unbiased regardless of the choice 
of vif(^). 

The C^^-type operations were performed using the 
conjugate-gradient algorithm, which was unprecondi- 
tioned for v and which used the preconditioner in Ap- 
pendix B of [13 with Zspiit = 64 for g. The traces to 
obtain Fa a are performed via the Z2 stochastic trace 
method [73]. 



C. Statistical errors 

Simulations are frequently used in CMB studies in or- 
der to (i) estimate error bars; (ii) verify that the analysis 
pipeline correctly reconstructs simulation inputs such as 
maps or power spectra; and (iii) determine the effects of 
possible systematic errors in the data. We use simula- 
tions for all three of these purposes, although only (i) is 
considered here (ii and iii are considered in Section [VI|1 . 
There are two major types of simulations, one with a 
circular beam and one with a toy model of an elliptical 
beam. The elliptical beam simulation is believed to be 
of higher fidelity, but we do a circular-beam analysis in 
order to understand whether beam ellipticity is an im- 
portant effect for this analysis. 

The circular-beam simulations described here are sim- 
ilar to those of Ref. 27], except that they now produce 
18 V and W band maps (one for each of the 3 years and 
6 DAs) instead of the 8 Q, V, and W band maps pro- 
duced in [13]. Additionally, the updated 3-year beam 
transfer functions B" of [351] are used in place of the 1- 
year values jTSj]. The simulation assumes a Gaussian 
unlensed CMB T and convergence k with power spectra 
generated by CMBFast [tJ] for the fiducial cosmology. 
(The convergence field we generate is Gaussian but based 
on the nonlinear power spectrum. However in our anal- 
ysis we will restrict to the linear regime, in which case 
Gaussianity should apply.) We also generate elliptical- 
beam simulations. The methodology for these simula- 
tions is described in Section IVC of Ref. [27*] and will 
not be repeated here; it consists essentially of taking the 
m = ±2 spherical harmonic components of the beam map 
for each DA and rotating them through a toy model of 
the WMAP scan strategy. We have of course updated the 
beam spherical harmonic coefficients based on the 3- year 
WMAP beam maps [33] ■ 
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V. RESULTS AND SIGNAL AMPLITUDE 

The cross-spectra between the large scale structure 
maps and the reconstructed convergence map are shown 
in Figure [31 The error bars are determined from 
cross-correlating the real LSS maps with 64 simulated 
elliptical-beam CMB lensing maps. 

Of greatest interest to us is the signal amplitude and its 
statistical significance. For each of the samples (LRGs, 
quasars, NVSS) we estimate the amplitude A defined as 
the ratio of the observed signal to the theoretical signal 
in the WMAP best-fit cosmology. That is, we fit 

Cf" (obs) =ACf"(th) (25) 

and determine the single parameter A. Since for each 
sample there are several bins in I, we require a weight 
(inverse-covariance) matrix in order to do a least-squares 
fit and find the best A. In principle the best covariance 
matrix for the Cf^ is that obtained from the simulations. 
However it is noisy and doing enough simulations to ef- 
fectively eliminate this noise would be computationally 
prohibitive. Therefore we have used the response matrix 
from the quadratic estimator procedure, Eq. ([20| . 
Since the matrix C*^^-* used in this procedure is not the 
true covariance matrix, the matrix used here is not 
necessarily the true covariance. For this reason it can 
only be used for weighting, and not for estimating the 
uncertainty in A. For the latter we must use the simula- 
tions. The amplitude is then determined from 

_^_ E^B[c^];iBgfr(th)gf;("bs) 

E^B[c.]ABCf;(th)c7f;(th) • 

Here the sums are over the Z-bins used in the fit, and 
{Ia, Ib} are the values of I at the centroids of the bin. In 
cases where we do not use all of the Z-bins, the weight 
matrices [C^,] are sub-blocks of corresponding to 
the /-bins used in the fit. The theory used is linear 
theory, and we set Zmax to correspond to 0.1/iMpc~^ at 
the lowest quartilc of the redshift distribution. In equa- 
tion form, Zmax = kmax/DA,25 where fcmax = 0.1/iMpc~^ 

and Da,25 is the comoving angular diameter distance to 
the 25th percentile of the bias-weighted redshift distri- 
bution f{z) (see Paper I). This maximum value of / is 
107 (LRGs), 269 (quasars), and 186 (NVSS). Because we 
exclude entire Z-bins if any range of multipoles is above 
the cutoff, this means we accept 5 bins (Zmax = 100) for 
the LRGs, 10 bins (/max = 250) for the quasars, and 8 
bins (/max — 175) for NVSS. If we combine the LRGs, 
quasars, and NVSS, the weighting procedure of Eq. 
gives relative weights of 0.239, 0.395, and 0.365 for the 
three samples, respectively. 

The results from this procedure are shown in Table [U 
the error bars are again determined by cross-correlating 
the real LSS maps with the 64 simulated CMB lensing 
maps. We have also shown in the table the results that 
are obtained by changing fcmax from its fiducial value of 
0.1 to 0.05 or 0.15 hMpc~^. As one can see these changes 



in the analysis result in < la changes in the amplitude 
A. 

The primary result of this paper is the standard pro- 
cedure using the TT bands (i.e. the signal averaged over 
the V and W bands). This is A = +1.06 ± 0.42, i.e. for- 
mally a 2.5(7 signal. Note that the relative weighting of 
LRGs, QSOs, and NVSS is determined from the response 
matrix, Eq. (|26p . If we used weights from the simula- 
tion covariance instead the amplitude would instead be 
-1-1.08 ± 0.42, a negligible change. 

We have also computed the mean values of A ob- 
tained from the simulated CMB lensing maps; these are 
-1-0.01 ± 0.10, +0.08 ± 0.09, and -f 0.03 ± 0.07. These are 
all consistent with zero at the la level. (We have quoted 
error bars on the mean of 64 simulations, which are a fac- 
tor of -\/64 = 8 smaller than the error bars on the data.) 
Since the LSS maps in our analysis are real instead of 
random, these mean values should be considered a test 
of whether the galaxy-convergence cross correlation that 
we observe is due to some feature of the galaxy map being 
aligned by chance with the CMB mask. 

A signal that is formally 2.5cr may not be statisti- 
cally significant if the noise distribution (value of A un- 
der the null hypothesis of no lensing) is non-Gaussian. 
We have therefore computed the skewness and kurto- 
sis of the Monte Carlo values of A. For the com- 
bined signal we find (A^) / {A'^f^'^ = +0.15 ± 0.48 and 
{A^)/{A^f = 2.46 + 0.61, which for a Gaussian distri- 
bution should be and 3 respectively. Thus with the 64 
simulations there is no evidence for any departure from 
Gaussian! ty. 



VI. SYSTEMATIC ERRORS 

This section considers possible systematic errors re- 
sulting from instrumental or algorithmic effects. Galac- 
tic emission, and some general tests on the reconstructed 
map. The extragalactic foregrounds require separate con- 
sideration, and discussion of them is deferred to Sec- 
tion |VlTl 



A. Calibration and beam ellipticity 

We have tested the calibration of the lensing estima- 
tor using the simulations from Section HV CI Aside from 
being an important test of the code, we note that the cal- 
ibration of the lensing estimator could differ from unity 
for several reasons. First, Eq. (I14p was derived assuming 
the small-deflection limit, i.e. it is constructed only to 
first order in the deflection angle d. The simulation does 
not use any Taylor expansion in d. Secondly, there are a 
number of places in the reconstruction and cross-power 
estimation codes where interpolations, /-cutoffs, and iter- 
ative C^^-type operations are used, and it is important 
to test whether these have introduced calibration biases. 
Finally, the beam ellipticity causes the effective transfer 
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TABLE L The amplitude A of the lensing cross-correlation signal normalized to that predicted for the WMAP cosmology 
(Eq. [25}. The first four rows show the "standard" fit as described in the text, following rows show the consequences of 
modifying this procedure. 



Fit type 


WMAP bands 


A (LRGs) 


A (QSOs) 


A (NVSS) 


A (combined) 


Standard 


TT 


+0.72 ± 0.76 


+1.20 + 0.73 


+1.11 + 0.52 


+ 1.06 + 0.42 


Standard 


VV 


+0.42 ±0.85 


+0.48 + 0.83 


+1.28 + 0.64 


+0.92 ± 0.48 


Standard 


VW 


+0.71 + 0.78 


+1.33 + 0.76 


+1.08 + 0.54 


+ 1.07 + 0.44 


Standard 


WW 


+1.16 + 1.19 


+1.80 + 0.91 


+0.98 + 0.66 


+ 1.21 + 0.54 


femax = 0.05/1 Mpc"^ 


TT 


+1.51 + 1.12 


+1.56 + 0.80 


+1.36 + 0.60 


+1.43 + 0.50 


fcmax = 0.15/1 MpC"^ 


TT 


+0.48 + 0.57 


+1.03 + 0.73 


+0.99 + 0.52 


+0.88 + 0.39 




LRGxWW 



QSOxWW 



10 100 



FIG. 3: The LSS-convergence cross-spectra. Each column represents a different large scale structure tracer (LRGs, quasars, 
and NVSS), and each row represents a different frequency combination (the first row, TT, is frequency-averaged; the other 
rows are VV, VW, and WW). The horizontal dashed line is zero, and the thick line shows the theoretical signal for the fiducial 
cosmology and the best-fit LRG bias/redshift distribution (see Paper I). The last data point in each plot is actually two /-bins 
that have been combined so that their error bars do not dominate the vertical scale of the plot. 
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TABLE II: The calibration of the lensing estimator, i.e. the 
mean value of Aout/Ain obtained from a set of 64 simulations 
with the fiducial cosmology, and the la error on this mean 
value. For perfect calibration this factor is 1. 



Beam 




^out/-^in 




type 


LRG 


QSO 


NVSS 


Circular 


0.930 ± 0.035 


0.979 ± 0.026 


1.059 ±0.019 


Elliptical 


0.953 ± 0.038 


0.959 ± 0.024 


1.040 ± 0.021 



function Bi to depend not just on the magnitude of the 
wavenumber I but also on its direction, which is not taken 
into account in constructing Eq. (jl4p . 

We have used simple simulations to estimate these ef- 
fects. In each case a Gaussian k map and corresponding 
lensed temperature maps were constructed, as described 
in Section IVC of Ref. [13] ■ These maps were beam- 
smoothed, and we added Gaussian uncorrelated noise 
(based on the WMAP A'obs values), which is accurate 
for the high multipoles where noise is significant. They 
were then fed through the reconstruction pipeline to de- 
termine v'^'^ . A galaxy field was generated with the cor- 
rect cross-spectrum Cf^ by setting gim = {Cf^ / Ci'^)Hiim- 
(In principle there is also an uncorrelated contribution to 
the galaxy field, i.e. we should be adding an additional 
contribution with power spectrum Cf^ — Cf'"^ /C^'^; how- 
ever when computing cross-correlations this term will av- 
erage to zero, and including it merely adds noise to the 
simulation results. We used the theoretical Cf^/C/"" for 
these simulations, based on the bias and redshift distribu- 
tions for the fiducial cosmology from Paper I.) The input 
power spectra are those from the fiducial cosmology, and 
the output spectra Ci'^ are fit to construct an amplitude 
A, just as with the real data. Since the input lensing 
amplitude is Ajn = 1, the output Aout can be used as 
a test of the calibration of the lensing pipeline. Ideally 
^outMin = 1- The results from this test are given in 
Table [IT] for both circular-beam and elliptical-beam sim- 
ulations; the elliptical-beam case is closer to reality, but 
the circular-beam case is useful for separating out the 
possible contributions to the calibration error. Note that 
the calibration bias ^out/^in — 1 is within 2a of zero in 
most cases, the exception being the case of the circular 
beam in cross-correlation with NVSS. The implied cali- 
bration error for this case is 6%, which is still much 
smaller than the statistical error. 



B. Power spectrum of reconstructed map 

One test of the reconstructed map is to compare its 
power spectrum to simulations. This test is useful be- 
cause it is sensitive to any source of spurious power in 
the lensing reconstruction such as foregrounds, unmod- 
eled beam ellipticity, or 1// noise. It is also impor- 
tant because the error bar on the LSS-convergence cross 
correlation depends on the convergence power spectrum 



(including noise and systematics), so it is essential to 
compare the measured convergence power spectrum to 
simulations in order to establish the reliability of the 
simulation-based error bars. 

In Figure IH we show the pseudo-C/ power spectrum 
of the reconstructed lensing field with the Kp052 mask. 
This power spectrum was derived via the formula 



pseudo-Cf"" = 



jv{h)-Yf:^{h)d^h 



{21 + l)i?f/sky 



(27) 



where are the vector spherical harmonics. As can 
be seen from Eq. (jl4p . this is simply the power spec- 
trum of the longitudinal component of v, divided by Rf 
to convert it to a convergence power spectrum. Note 
that since our objective here is simply to compare the 
amount of power in the data with the simulations, we 
have not done any noise subtraction, and we have di- 
vided by /sky instead of doing a full deconvolution of the 
Kp052 mask. The pseudo-C;s in the figure have been 
binned in spacings = 10. We have also shown this 
power spectrum for the individual- frequency maps v^^, 
v^'^, and v'^'^; since the power spectrum of the con- 
vergence is noise-dominated, it is of course lowest in the 
frequency-averaged map v. 

For comparison, we repeated the same analysis for the 
Kp05 cut, i.e. without masking the point sources. This 
resulted in excess power that is not present in the simu- 
lations; this is especially the case in the map, which 
has the highest point source contamination. For this rea- 
son the Kp05 mask was not used in the rest of our anal- 
ysis. 



C. 90° rotation test 

An important systematics test in lensing analyses us- 
ing galaxies as sources is to rotate each galaxy ellipticity 
by 45 degrees and look for a lensing signal. This is com- 
monly known as the "_B-mode test" since the rotation by 
45 degrees interconverts the i?-mode shear pattern pro- 
duced by lensing and the S-mode pattern that should 
be zero in the absence of systematics. For our analy- 
sis, the reconstructed CMB lensing field is a vector with 
spin 1, hence the analogous test here is to rotate v by 
90 degrees and re-measure Cf^. This is done in Figure [5] 
for each of our samples. The for ^ signal using 
the covariance matrix from rotated simulations is 7.85 
for the LRGs, 14.55 for the quasars, and 17.32 for NVSS 
for 13 degrees of freedom. One can also fit for the am- 
plitude Aqq of these rotated cross-power spectra using 
the theoretical lensing signal as a template and the same 
range of I as for the main fit [analogous to Eq. (pS)) . ex- 
cept that this time one expects to get zero]. This gives 
Ago = -0.61 ± 0.72 for the LRGs; Ago = +0.77 ± 0.62 
for the quasars; and = -t-0.16 ± 0.33 for NVSS, with 
the error bars determined from the rotated simulations. 
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Power spectrum of reconstructed lensing field with Kp052 mask 
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FIG. 4: The pseudo-C; power spectrum of the reconstructed 
lensing map (Ea. l27p for the Kp052 cut. These power spectra 
are dominated by noise at aU scales. The top panel shows the 
power spectrum of the full reconstruction, whereas the other 
panels show the power spectra of the v^^, v^^, and v'^^ 
lensing maps. The solid line shows the data, and the dotted 
lines show five simulations. Note that with the Kp052 sky cut 
the simulations reproduce the power spectrum of the data to 
within several percent. Similar plots for the Kp05 cut show 
significant excess power. 



These are indeed consistent with zero (the largest value 
is for the quasars which are 1.24cr from zero). 



D. Frequency dependence of signal 

One simple but important test for the CMB lensing 
origin of the signal we observe is the frequency depen- 
dence. If CMB lensing is the correct explanation then 
the signals in the VV, VW, and WW lensing maps will 
be identical (aside from instrument noise). If instead the 
signal is due to some other effect, such as point sources 
correlated with the LSS tracer, then a different frequency 
dependence is expected. In general for each spectrum we 
can construct the ratio of thermodynamic temperatures 
Rv/w = T{V)/T{W). The lensing amplitudes A ob- 
tained for each of the lensing maps should then be in 



the ratio R 



v/w 



:1:R 



V/W 



for the VV:VWWW maps. 



respectively. The most worrying foregrounds are flat- 
spectrum radio sources (Ry/w = 2.1 for a = 0), infrared 
sources {Ry/w = 0-46 for a = 3.5), and the thermal 
Sunyaev-Zel'dovich (tSZ) effect {Ry/w = 1.16). By com- 
parison, frequency-independent CMB fluctuations should 



have R 



v/w 



1. 



The simplest way to set a constraint on Ry/w is as 
follows. For each LSS sample, we obtained in Section IVl 
the amplitudes Ayy, Ay\y, and Awwi which form a 3- 
component vector A, and their 3x3 covariance matrix 
Cova is known from the Monte Carlo simulations. We 
also know that given Ry/w the vector A should be paral- 



lel to the template vector t{Ry/y/) 
Then we construct the x^: 



{Ry/W^ l,i?y/^). 



X — min[A 

V 



yt{Ry,w)V[CoYA\-\^ - y^{Ry/w)] 



(28) 

The curves for each sample (LRG, quasar, and NVSS) 
and the combined case are shown in Figure [6l For the 
combined case, we find Ry/w = 0.87tJ]:^:^ (2cr, deter- 
mined by setting x^ - Xmin = 4)- 

For each of the LSS samples used in this paper, the fre- 
quency dependence of the LSS-convergence correlation is 
consistent with a blackbody spectrum {Ry/w = !)■ In 
particular for the combined LSS sample, the frequency 
dependence disfavors the possibility that the signal is 
due to flat-spectrum radio sources at 2.3cr. With the fre- 
quency dependence tests alone, we cannot rule out the 
possibility that the signal is due to infrared sources or 
tSZ haloes; our constraints on these contaminants are 
discussed in Section IVlIl 



E. Galactic foreground 

Galactic foreground emission is an important potential 
source of systematic error for CMB experiments, includ- 
ing those aiming to study lensing. It arises locally and 
thus cannot correlate with LSS, but it can correlate with 
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LRGxTT(rot) QSOxTT(rot) NVSSxTT(rot) 
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I t I 

FIG. 5: The LSS-convergence cross spectrum Cf* obtained by rotating v 90 degrees before computing the cross-correlation. 
This should be zero in the absence of systematics. The three panels show the LRGs, quasars, and the NVSS sample. The solid 
lines show the signal expected for the unrotated maps in the WMAP cosmology. The last data point in each plot is actually 
two /-bins that have been combined so that their error bars do not dominate the vertical scale of the plot. 



Constraints on frequency dependence of "lensing" signal 
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Rv/W=T(V)/T(W) 

FIG. 6: Constraints on the frequency dependence of the lens- 
ing signal. The horizontal axis is the ratio of thermodynamic 
temperatures in the V and W bands, which is 1 for the CMB. 
The vertical axis is the obtained by fitting the values of 
the lensing amplitude A in the VV, VW, and WW frequency 
combinations. 



systematic errors in the LRG and quasar maps, such as 
stellar contamination or errors in the extinction correc- 
tion. Galactic effects can be treated in basically three 
ways: (i) one could search for frequency dependence in 
the maps; or (ii) one could feed an estimate of the Galac- 
tic foreground into the lensing estimator and correlate 
this with the real LSS maps. Of these, (i) is in prin- 
ciple the best if a wide range of frequencies, e.g. the 
whole WMAP range 23-94 GHz, is available. Unfortu- 
nately the lower-frequency bands do not resolve the high 
multipoles so this approach has limited applicability, and 
in any case the error bars associated with this approach 
are large. Therefore we also use (ii), which gives much 
tighter constraints. 

A model for the Galactic foreground must include sev- 
eral known microwave-emitting components. In the V 



and W bands, thermal emission from dust and free-free 
radiation from electron-ion collisions are the dominant 
foregrounds. At the lower WMAP frequencies there is 
also a soft synchrotron component. In addition, there is 
an "anomalous" component that dominates at the lower 
WMAP frequencies and appears to be significant in V 
band, which may be due to rotational or magnetic dipole 
emission from dust. 

We use a model for the Galactic foreground based on 
the analysis of [s^l ; this is the same model as was used in 
the first-year WMAP lensing analysis [27| . This model is 
composed of the "Model 8" [Sll, [s^l prediction for ther- 
mal dust emission, added to the prediction for free-free 
radiation based on Ha line emission [83| with the conver- 
sion factor of [s^l . The anomalous microwave emission is 
included by multiplying the thermal dust template by a 
factor proportional to T^^^^, which provides a good phe- 
nomenological fit [83|. The soft synchrotron component 
is negligible in V and W bands, so we do not include 
it. (In any event, the usual soft synchrotron templates 
used in CMB work, namely the 408 MHz Haslam maps 
(ssl [soj , have a ~ 50 arcmin beam that filters out the 
high multipoles used in the lensing analysis. Use of this 
map as a foreground indicator for lensing or any other 
high-/ application would give misleading results.) 

Since the Galactic emission is not correlated with the 
CMB, we may estimate the Galactic contamination to 
the lensing amplitude AA by feeding the Galactic fore- 
ground model maps through the lensing pipeline in place 
of the actual WMAP maps. We computed this for each 
of the frequency combinations (VV, VW, and WW) and 
each sample (LRG, quasar, and NVSS); the worst value 
of the contamination \AA\ was 9 x 10^'*. This is negli- 
gible compared to our error bars and hence we make no 
attempt to correct for it. 



VII. EXTRAGALACTIC FOREGROUNDS 

This section considers the major extragalactic fore- 
grounds that can contaminate the lensing signal: point 
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sources; thermal SZ fluctuations; kinetic SZ; and the 
Rees-Sciama efi^ect. The point sources can be subdi- 
vided into "radio" sources, which emit via synchrotron 
and free-free radiation and are brightest (relative to 
the CMB) at low frequency, and "infrared" sources, i.e. 
dusty galaxies. We rely mostly on a combination of fre- 
quency and spatial dependence of the foregrounds to sep- 
arate them from lensing. The brightest radio sources 
are usually active galactic nuclei (AGN), although star- 
forming galaxies also emit synchrotron and free-free. 
Note that it is possible for the same object to be both a 
radio and an infrared source. 

This section is organized as follows. The spatial depen- 
dence of foregrounds is discussed in Section IVII Al where 
we show that the bias in the lensing signal depends on the 
galaxy-foreground-foreground bispectrum. This bispec- 
trum is decomposed using the halo model and we show 
how to constrain each term using the galaxy-temperature 
cross-spectrum and physical constraints on the point 
sources such as flux cuts. Most of the rest of the section is 
based on this framework. We apply the constraints from 
Section rVlI Al to radio point sources in Section lVlI Bl and 
infrared sources in Section IVII CI We then investigate 
the tSZ effect (Section IVIID[) . Several small effects are 
considered in the appendices: Appendix [Cl considers the 
kSZ and ISW/Rees-Sciama effects, and Appendix |D] con- 
siders correlations between different foreground compo- 
nents. The arguments used to constrain the extragalactic 
foregrounds are summarized in Section IVII El where we 
discuss what we believe are the most significant strengths 
and weaknesses and what could be done to improve them 
in the future. 



A. Point sources: spatial correlation-based tests 

Point sources are the most serious foreground for the 
WMAP lensing analysis. They introduce two major con- 
cerns: one is that the artifacts they produce in the lensing 
map can correlate with large-scale structure, thereby in- 
troducing a spurious cross-correlation signal; and they 
can add power to the lensing map and thus increase the 
error bars on the cross correlation Cf^. Neither effect 
is taken into account in the simulations. However the 
comparison of the lensing power spectra in Section [VI Bl 
demonstrates that the point sources do not add signifi- 
cantly to the power in the v map and hence do not in- 
crease the errors f7(Cf^). 

We have several methods of constraining the bias in 
C^'" due to point sources: (i) tests based on the fre- 
quency dependence of the signal; (ii) tests based on 
spatial correlations; and (iii) tests based on external 
foreground maps. Of these, (i) was considered in Sec- 
tion lVI D1 and (iii) will be discussed using the IRAS maps 
of dusty sources in Section lVlI C 21 This section will focus 
on the methodology for (ii) , the tests using spatial corre- 
lations. This methodology is applied to radio sources in 
Section IVlI Bl and infrared sources in Section IVlI CI 



Since our estimates of Cf^ are obtained by cross- 
correlating the LSS map with a "convergence" map 
that is quadratic in temperature, it follows that any 
point source contamination must enter via the cross- 
bispectrum B^^^Ji that contains the LSS tracer and two 
temperatures. The amount of contamination to the Cf" 
can be written as 

Acr - i?r' E ^u'i"Bfi^;j;{ps); (29) 
I'l" 

see Appendix [X] for a derivation and Eq. (jA2p for the 
definition of the coefficient J-wi". 

We will describe the point source contributions to the 
cross-bispectrum Bf^f^Ji in the context of a halo model. 
Note that since most haloes are unresolved by WMAP, 
the haloes are described entirely by the number of "galax- 
ies" (LRGS, quasars, NVSS sources) they contain, their 
microwave flux in each band (in thermodynamic fiK sr), 
and their 2- and 3-point correlation functions. In most 
cases we ignore that the haloes are resolved since only 
the nearest haloes are resolved and in any case we have 
directly tested that an NEW profile leads to a lower 
contamination to the lensing signal than a delta func- 
tion with the same integrated flux. Since very little in- 
formation is available on the halo occupation distribu- 
tions (HODs) for quasars, NVSS sources, and microwave 
sources, our analysis will be based to the maximum ex- 
tent possible on general properties of the halo model; in 
particular we will avoid parameterized HODs entirely. 

Within the context of the halo model, the point source 
contributions to the cross-bispectrum Bf^^^J^ can have four 
contributions: 

1. The "1-halo" term, which arises when the two fac- 
tors of temperature come from point sources in the 
same halo as the galaxy. 

2. The "2a-halo" term, in which the two factors of 
temperature come from point sources in the same 
halo, but the galaxy is in a different halo. 

3. The "2b-halo" term, in which there is one point 
source in the same halo as the galaxy and the other 
point source lives in a different halo. 

4. The "3-halo" term, in which the galaxy and both 
point sources live in three distinct halos. 

Note that there can be no galaxy-point source-primary 
CMB cross-bispectrum since the primary CMB resides at 
z ~ 1100 whereas the galaxies and point sources are at 
much lower z. 

The formulas for Bf^j^J] are simply stated here so that 
we can concentrate on the physical assumptions required 
by our argument; derivations for the 1- and 2-halo terms 
can be found in Appendix [B] while the 3-halo term is de- 
rived by the usual Limber argument (see e.g. Ref. [87|). 
The number of galaxies N and microwave flux of a 
halo are of course dependent not only on the intrinsic 
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properties of the halo but also on the redshift, since flux 
declines as the luminosity distance squared (times the k- 
correction) and the selection criteria for our samples are 
redshift-dependent. For brevity we will not write this de- 
pendence explicitly. The 1- and 2-halo terms make use 
of the maximum unmasked point source flux i^max and 
the maximum number of galaxies per halo Ny^nx- We 
use the notation p{F) dF is the comoving 3-dimensional 
number density of haloes emitting microwave flux be- 
tween F and F + dF and p{N) to denote the comoving 
3-dimensional number density of haloes with A'^ galax- 
ies. The 2-dimensional (projected) number densities will 
be denoted n2D{F) and n2D{N). We will also use the 
notation 



fix) 



d7 

Y^Np{N)b{N)^—J{z 

N ^ 



(30) 



for the bias-weighted comoving distance distribution of 
the galaxies. 

We will consider the 1- and 2a-halo terms together, 
since they are closely related, and then the 2b- and 3- 
halo terms separately. We will use two separate ways 
of estimating the 1- and 2a-halo terms. Method I uses 
the cross-power spectra Cf"^ obtained in Paper I, and 
crudely speaking the basic idea is that for 1- and 2a- 
halo terms, the conversion factor from power spectrum 
to bispectrum cannot be greater than the maximum flux 
^max- Method II looks at the bispectrum of the galax- 
ies and two frequency-differenced maps to constrain the 
F^-weighted correlation function of the galaxies and the 
point sources. Generally for types of sources that are very 
faint, i.e. where -Fmax is below the WMAP noise level per 
beam. Method I gives the tighter constraint, whereas for 
radio sources (which can be extremely bright) Method II 
is superior. 



1. 1- and 2a-halo terms: General considerations 

The 1-halo contribution to the cross-bispectrum is sim- 
ply the Poisson bispectrum, 

/N 
dFj2^nMN,F), (31) 
N ^3 

where F is the halo flux, N is the number of galaxies in 
the halo, fig is the mean number of galaxies per steradian, 
and n2D{N, F) dF is the number of haloes per steradian 
containing N galaxies and with flux between F and F + 
dF. The coefficient rjunn is given by 



(2^ + i)(2/- + i)(2/" + i) / I I' r \ 
m'l" - \i ^ 0' ^^'^ 



and simply encodes the products of spherical harmon- 
ics. Usually Eq. (pij) would be written in terms of 3- 
dimensional densities and then integrated over redshift 



(or comoving distance x) using the Limber approxima- 
tion. This is useful for theoretical predictions since the- 
ory usually gives 3-dimensional densities and biases of 
haloes as a function of mass and redshift, and does not di- 
rectly give results in projection. However in our case both 
the desired constraint (on B^^^^) and the given informa- 
tion {C^^ measured from Paper I) are already projected 
onto the sky, so we can derive our constraints entirely 
with Eq. pip without considering the 3-dimensional dis- 
tribution of the point sources. 

The analogous equation to Eq. (pij) for the 2a-halo 
term is 



<"(2ah) 



" AT 



NFF' 



mi'i' 

xn2D{N)n2D{F)Ci{N;F), 



(33) 



where n2D{N) is the 2-dimensional density (objects per 
steradian) of haloes with N galaxies; n2D (F) dF is the 
2-dimensional density of haloes emitting flux between F 
and F + dF; and Ci{N;F) is the cross-power spectrum 
between haloes containing N galaxies and haloes with 
flux F. This cross-power does not include shot noise. 
Once again for the sake of ab initio theoretical predictions 
it would be most useful to write this as a Limber inte- 
gral, containing 3-dimensional instead of 2-dimensional 
densities and involving the 3-dimensional cross-spectrum 
between different types of haloes; but for our purposes 
we can work entirely in the 2-dimensional space. 

This implies a combined l-|-2a-halo contribution to the 
cross-bispectrum. 



i?f,?;r(l + 2ah) = m'l" fdFj2 



NF' 



[n2D{N, F) 



+n2D{N)n2DiF)CiiN;F)]. (34) 

The corresponding bias in the cross-power is obtained 
from Eq. 



ACr(l + 2ah) = [Y.^^ 

AT 



H'Wii'i" 



[712d{N,F) 



N -a 

+n2D{N)n2D{F)Ci{N;F)]. (35) 

Note that /' and /" enter only in the summation in 
parentheses, a fact that will prove to be extremely valu- 
able. This "compartmentalization" is a direct result of 
the pointlike nature of the sources under consideration. 
It is convenient to deflne rps{l) to be the combination 
in parentheses (this is consistent with the definition in 
Ref. [131), so that 



ACT(l + 2ah) = 



Ri 



NF' 



[n2DiN, F) 



N -"9 

+n2D{N)n2D{F)Ci{N;F)]. (36) 
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The point source responsivity rps{l)/Ri is negative for 
I < 269. This is a measure of the spurious lensing signal 
that is produced by a point source as a function of scale. 
Note in particular that if the galaxies and point sources 
are positively correlated, i.e. 

J2N[n2D{N,F)+n2DiN)n2D{F)Ci{N;F)]>0, (37) 

N 

then the bias AC;^'^(l+2ah) must have the same sign as 
rps{l)/Ri- Note that for / < 269 we have rps{l)/Ri; this 
is simply a statement that the long-wavelength {I < 269) 
modes of the reconstructed lensing map have a negative- 
convergence artifact around each point source. Since our 
fits to the lensing amplitude A use only the / < 269 mul- 
tipoles, it follows that it is not possible for the observed 
lensing signal to be due to the 1- or 2a-halo point source 
terms (i.e. artifacts in the v map that arise from indi- 
vidual sources rather than interference between sources 
clustered with each other). It is possible for the ampli- 
tude A to be biased downward (at some level it must be) 
and this must be considered in any cosmological analysis. 



2. 1- and 2a-halo terms: Method I 

Our objective here is to set a constraint on the possible 
contamination, Eq. To do this, we consider the 

cross-spectrum of the galaxies and the temperature, Cf^, 
coming from the point sources. This is a relatively simple 
calculation that has only 1- and 2-halo contributions; the 
answer is 

Cf (ps) = f dFY,^[n2DiN,F) 

+n2D{N)n2D{F)Ci{N-F)]. (38) 



3. 1- and 2a-halo terms: Method II 

We now consider another way to constrain the 1- and 
2a-halo terms. This time, instead of using a maximum 
flux Finaxj we attempt to directly measure the integral 
in Eq. ([36)) using a quantity quadratic in the CMB tem- 
perature and linear in the galaxy density. We do this us- 
ing frequency-differenced maps, which renders our point 
source analysis completely insensitive to biases produced 
by lensing of the CMB, and to noise produced by the 
primary CMB anisotropy. This method is useful only for 
the radio point sources since for these Method I is not 
very constraining (due to large Fmax)- 

The basic idea is that we produce two frequency- 
difference maps andD^^). In the simplest version of 
the method, the map is the difference of the WMAP 
temperature maps T{Ka) — T(y), both smoothed to the 
resolution of the Ka-band map. The map Z?^^^ is simi- 
larly constructed from T{Q) — T{W). We will denote the 

effective (post-smoothing) beam transfer functions B^^' 

(2) 

and _Bj . We then multiply the two maps to produce 
D{n) = L»(i)(n)i:)(2)(n). The idea is that each point 
source produces a feature in -D(n) whose intensity is pro- 
portional to F^ and various factors that involve the fre- 
quency dependence. Specifically, if the input map con- 
tains a single point source, then the output map intensity 
can be shown to equal 



(40) 



where fii is the position of the point source, Fy^i is its 
V-band flux (in units of /iKsr), r/ is the transfer function 



I'l" 



(2r + l)(2r + l) ll'l"] ^ (41) 

4^ ' ' I / ' ^ ^ 



We now note that the integral over dF in Eq. f36|) 
is identical to Cf^(ps) except for an additional factor 
of F. If we make the assumption that the fluxes are 
less than some i^max, and also assume that the inte- 
grand is nonnegative (Eq. [57|) . it follows that the in- 
tegral over dF in Eq. p6|l is bounded between and 
^maxC'f"^(ps). Then we can conclude that (i) the con- 
tamination ACf'^{l + 2ah) has the same sign as rps{l), 
and (ii) the contamination is bounded by 



iAcr(i 



2ah)| < 



Ri 



Cf (ps). 



(39) 



This simple inequality, which allows us to estimate the 
point source contamination, relies only on the existence 
of a bound on the point source flux, and the require- 
ment of Eq. (|37p . i.e. that galaxies and point sources are 
positively correlated (where the correlation includes the 
Poisson term) . Note that aside from the positivity of the 
correlation no assumption is made about the halo mass 
function, bias, environment dependence or lack thereof, 
etc., so we believe this represents a very robust bound. 



and the spectral-dependent terms are if'^^ = {Fxa^i — 
Fv.i)/Fv,i and similarly for K'-^\ (This formula follows 
immediately from the multiplication in real space of the 
two beams.) If the map contains multiple point sources 
then in addition to terms of the form Eq. (|40p there are 
interference terms involving combinations of two point 
sources (with different sources contributing to D^^^ and 
\Ye may then do a cross-correlation with the 
galaxies, Cf^. The l-|-2a halo contributions (i.e. leaving 
out the interference contribution to D) to this are 

Cf ^ (1 + 2ah) = K^^^Ti J dF n2D [F)F^Ci {g,F), 

(42) 

where Ci(g,F) is the angular cross-spectrum between 
galaxies and haloes of flux F. This can be expanded 
as 

/Af p2 
dFY,-^[ri2D{N,F) 
N ^9 

+n2D{N)n2D{F)Ci{N;F)], (43) 
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so that 



ACT(l + 2ah) 



2ah) 



Ri KWK(^)ti 



(44) 



We use this equation to estimate the l+2a halo point 
source contamination of Cf^. 

In this method it is important to remove the Galac- 
tic foreground, particularly since one is working with 
the low- frequency bands where the Galaxy is brightest. 
There are two ways to do this. One is to subtract the 
Galactic foreground model (Section IVI Ep . Alternatively, 
instead of directly multiplying the maps D^'^\ one can 
high-pass filter them, which should remove more of the 
slowly varying Galactic signal. The latter method also 
has the advantage that the temperature offsets in the 
different bands is not needed (WMAP does not measure 
these directly and used a simple model of the Galaxy 
to set these offsets). A simple, local high-pass filter is 
to subtract the map smoothed by a Gaussian of width 
CTi = 45 arcmin; this is equivalent to modifying the effec- 
tive beam in accordance with 



B 



(45) 



Thus Ti is different for the two methods. 



^. 2b-halo term 



The analogous equation to Eq. pTjl for the 2b-halo 
term is 



Bf^^,{2W) = m'l" I dF I dF'J2 

J J 



NFF' 



xn2D{N, F)n2D{F') 

x[Cv{N,F-F') + Ci.{N,F-F% (46) 

where n2D{N) is the 2-dimensional density (objects per 
steradian) of haloes with N galaxies; n2D{F)dF is the 
2-dimensional density of haloes emitting flux between F 
and F+dF; and Ci{N, F; F') is the cross-power spectrum 
between haloes containing N galaxies with flux F and 
haloes with flux F' . The spurious contribution to the 
galaxy-convergence correlation is obtained from Eq. ((29)) : 



ACr(2bh) = 2R-^Y.^ii'i"'nii'i- 
I'l" 



fdF fdF'Y: 



NFF' 



N "9 

xn2D{N, F)n2D{F')Cv (N, F; F'). (47) 

Here we have removed the Ci"{N, F; F') term and dou- 
bled the result because of the symmetry of J-' and rj under 
interchange of I' and I". 

We would like to place a bound on AC;^'"(2bh) us- 
ing the galaxy-temperature cross power spectrum, just as 
was done for the 1- and 2a-halo terms in Method I. There 



are two simple ways of doing this: one based on the max- 
imum flux i^maxj and the other based on the maximum 
number of galaxies per halo Nmax- The former is model- 
independent and is preferred, except for the radio sources 
for which it docs not give an interesting constraint due 
to large Fmax. 

The -Fmax method works by using the triangle inequal- 
ity in Eq. (|T7|) . and then noting that F < i^max: 



iAcr(2bh)i < 2i?r^E 



vqii'i'- 



N 



NF F' 



xn2D{N)n2D{F')Ci,{N-F'). (48) 

The last two lines in this equation are - aside from the 
factor of i^max ~ equal to the 2-halo galaxy-temperature 
spectrum. This must be less than or equal to the actual 
galaxy-temperature spectrum since the 1-halo term must 
be nonnegative, so we may write 



iAcr(2bh)i < 2i?r'E 



xF„,axCf (ps) 



(49) 



The last cross-spectrum was directly measured out to 
I' = 600. It is extrapolated to higher I' assuming 
Cf,^ ^constant, which is conservative since the 2-halo 
Ci' decreases at high multipole [due to the declining P{k) 
over the relevant range of scales]. 

The model-dependent iVmax method works as follows. 
Written in terms of angular correlations, the 2-halo con- 
tribution to the point source power spectrum is 

Cr^(ps,2h) ^ J dF j dF' FF'n2D{F) 

xn2D{.F')Cv{,F-F'). (50) 

We then have 



iAcr(2bh)i < 2i?r'E 



E^' 



u'i"Vu'i" 



max ^TT 



(Ps,2h), 



(51) 



where A'max is the maximum number of galaxies found 
in a halo and we have assumed the positivity condition 



J2NCiiN,F;F')>0. 



(52) 



N 



On the other hand, written in terms of a Limber integral, 
Cf^(ps,2h) is 

Cr(ps,2h)= / ^[bTf'{x)fPUk) (53) 
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where Piin(fc) is the Unear power spectrum at fc = 
l/2)/r-, 



T'ix) ^ 



dFFp{F) 



(54) 



is the mean contribution to the brightness temperature 
per unit comoving distance, and 



br = 



_ J dFFp{F)b{F) 
J dFFp{F) 



(55) 



is the effective (flux-weighted) bias of the microwave 
emitters. We note that T'{x) does not diverge at small 
distances even though the fluxes from certain haloes can 
become large because of the factor; this can be seen 
explicitly from the relation 



2kv 



r(l + ^)j 



com 



(56) 



where v is the observed frequency, j™™ is the comoving 
emissivity (in e.g. ergcm~^s~^ Hz~^ sr~^, where the cm 
are comoving and all other units are physical), and is 
Boltzmann's constant. 

This equation is not quite in a usable form, since the 
intensity and redshift distribution (T") and bias (&t) of 
the non-primordial background are still unknown. There- 
fore we have tried a different approach. Using the Limber 
approximation, and assuming linear biasing, the second 
(2-halo) term of Eq. ((38| can be rewritten as 

Cf > Cf^ (2h) 

= fdxr' [ dFY,^PiF)p{N) 



xb{F)b{N)PUk), 
= f ^f{x)bTf'{x)PUk). 



(57) 



Physically we expect the bias bx to vary slowly with red- 
shift. The same is also true of ji, and hence T'{x) (see 
Eq. [55)1 . If over the range of redshifts of interest it is a 
power law, bxT' {x) ^ (1 + then 



6tT'(x)U=o < 



XT 



Jdxr-^f{x){l + z)iPUko)' 



(58) 



where equality would hold if the right-hand side included 
only the two-halo term, Iq is any arbitrary multipole, 
and fco = {Iq + l/2)/r. In practice to increase signal- 
to-noise ratio, we average the right-hand side over all 
the multipole bins between 25 and /max (the multipole 
corresponding to fc = 0.1/iMpc^^ at the 20th percentile 
in distance). We construct bounds on the 2b-halo point 
source contamination by considering a range of exponents 
(. Of course, all of these integrals to the redshift range 
where there are galaxies; we use Zmax = 0.8 for the LRGs, 
2.7 for the quasars, and 3 for NVSS. 



For the quasars and the NVSS galaxies the constraint 
on the maximum number of galaxies per halo will be 
redshift-dependent. In this case, it is straightforward to 
show that Eq. ([51]) undergoes the replacement 



iVmaxQ ' (pS, 2h) 



dx 



NmUx)[bTT'{x)VPun{k), 



(59) 

i.e. iVmax is pulled inside the integral over comoving dis- 
tance (or redshift). 



5. 3-halo term 

The 3-halo contribution to the point source bispec- 
trum is different from the 1- and 2-halo terms in that 
we have found no simple model-independent argument 
to constrain it using the galaxy-temperature correlation. 
Instead we will have to use a more model-dependent 
method. The Limber equation allows us to express the 
galaxy-temperature-temperature angular bispectrum in 
terms of the 3-dimensional bispectrum. 



Bwi^{3h) = rjwi" I dxr 



dF j dF'Y,p{N)p{F) 

'' N 



NFF' 

xp{F') B{k,k',k"\N;F;F'), (60) 

nx 

where x is the radial comoving distance; r is the comov- 
ing angular diameter distance (equal to x for a flat uni- 
verse); p{N) is the comoving number density of haloes 
containing N galaxies; p{F)dF is the comoving num- 
ber density of haloes with microwave flux between F 
and F + dF; the wavenumbers are k = {I + l/2)/r, 
etc.; and B{k, k' , k"\N; F; F') is the 3-dimensional cross- 
bispectrum of the haloes with N galaxies, the haloes with 
flux F, and the haloes with flux F' . The problem is that 
the key quantities p{F) and B(k, k' , k"\N; F; F') are not 
known. 

To circumvent this problem, we will consider the linear 
bias prediction for the 3-halo term. That is, we set 

B{k, k', k"\N; F; F') = b{N)b{F)b{F')B2p{k, fc', k"), 

(61) 

where b represents the bias of the haloes and 
B2p{k,k',k") is the second-order perturbation theory 
bispectrum. This allows us to factor Eq. ([SO]) , yielding 



^ii'i" 



dx 



fix)[bTf'ix)?B2Pik,k',k"). (62) 



Using Eq. ([29]) . we may then write 



iAcr(3h)i < i?r'Ek"' 



, dx 



I'l" 



xf{x)[bTf'{x)?B2p{k,k',k") .{63) 

The factor bTT'{x) comes from the bound, Eq. just 
as for the 2b- halo term ( A^max method) . Like the 2b- halo 
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term, it is model-dependent. Fortunately, the bound is 
also very small so deviations of halo clustering from the 
simple perturbation theory predictions or more compli- 
cated redshift evolution of bTT'{x) will have little effect. 



B. Estimates of contamination: radio sources 

We now construct radio point source contamination 
estimates using the methodology of the previous section. 
For the 2b- and 3-halo terms we need to know the contri- 
bution to the cross-spectra Cf^ from radio sources in or- 
der to get bxf'. We will use only Method II for the l+2a- 
halo terms because for the bright radio sources Method 
I does not give an interesting constraint. 

The Cf ^ cross-powers for the LRGs, the SDSS quasars, 
and the NVSS sources were obtained in Paper I in the 
range up to Z = 600. Unfortunately they are noisy, with 
the noise dominated by primary CMB fluctuations, and 
on large angular scales the ISW effect dominates the 
cross-power. For this reason, we have used the differ- 
ences Cf'^(O) - Cf'^{V) and Cj>^{V) ~ Cf^(PF) as the 
major sources of information on point sources. We con- 
struct our estimates for the V-band lensing estimator, 
which will be the most contaminated by radio sources. 
For flat spectrum sources and the relative weighting of 
bands in Eq. (fT7|) . the ratio of contamination is 



ACf^iTT) 



= 0.56. 



(64) 



We take as our estimate of the radio point source cross- 
spectrum 

Cf (l/;rps) EE ai[Cf (Q)-Cf (F)] 

+a2[Cf (V^)-Cf (W^)], (65) 

with ai = 1.128 and 02 = —0.466. The ratio 02/01 — 
—0.413 is chosen to cancel any tSZ contribution, and the 
normalization was chosen to be correct for flat-spectrum 
radio sources with ^constant. If there are also steep- 
spectrum sources this is an overestimate of their effect: 
for oc 1^-°-^ we have C'f^(F; rps) = 1.66Cf'^(F; rps). 
Equation (|65p has some sensitivity to infrared sources: 
for F^ oc 1^3-5 we have Cf^{V] rps) = 0.008Cf^(F; dust), 
and the coefficient 0.008 rises to 0.146 if F^ (x v'^. Since 
this coefficient is positive we will overestimate the con- 
tribution of radio sources if the infrared sources are also 
significant. This cross-spectrum is obtained by construct- 
ing a WMAP map of the frequency combination 1.128Q- 
1.594V-I-0.466W and cross-correlating it with each galaxy 
map. Also in order to speed up the calculation we did 
not implement any weighting in the Cf'^ estimation, i.e. 
instead of full weighting we only project out the 

monopole and dipole from each map (and in the case of 
NVSS, the declination rings and Haslam map). We have 
computed error bars for each case by cross-correlating 20 
simulated LSS maps with the actual WMAP difference 



map (on large angular scales in WMAP the difference 
maps contain 1// noise and Galactic foregrounds, and 
hence are more difficult to reliably simulate than the LSS 
maps). The simulated LSS maps are Gaussian and have 
the correct autopower spectra as determined from Paper 
I. The Gaussianity of the LSS does not matter for get- 
ting the variance of Cf{V] rps), since the Cf^(F; rps) 
estimator is a linear function of the data. 

The 2b- and 3-halo terms can be constrained if we 
know the maximum number of galaxies per halo A^max 
(needed for the 2b-halo term only) and the redshift evolu- 
tion exponent C, = dln[&Tr'(x)]/dln(l -I- z). We estimate 
A^max for NVSS by examining each object and finding the 
number of neighbors (including the object itself) within 
6 arcminutes (~ lh~^ Mpc at z = 0.2). For NVSS there 
are no objects with > 12 neighbors, 1 with 12 neigh- 
bors, and 2 with 11 neighbors, so we set A^max = 12. 
This may be a large overestimate: in particular Poisson 
statistics predicts 0.1 objects with 12 neighbors and 0.9 
with 11 neighbors, so it is not clear that these groupings 
should be taken seriously. Nevertheless for our purposes 
all we need is an upper limit, so we use 12. For the 
quasars the maximum number of neighbors is 10, so we 
set A'lnax = 10. Again this should not be taken espe- 
cially seriously (for Poisson statistics there should be on 
average 1.2 objects with 10 neighbors). For the LRGs 
a similar argument gives A^max = 36. For quasars and 
NVSS the redshift distribution extends to higher z where 
we can improve the constraint on A^max- For example at 
z = 0.8 a l/i"-'^ Mpc radius corresponds to 2 arcminutes; 
thus for z > 0.8 we have replaced the above A'max values 
with the maximum number of neighbors within 2 arcmin. 
This is 5 for the quasars and 5 for NVSS. For NVSS we 
have also set A"max = 3 for z > 3, which is what we find 
with a 45 arcsec radius {\h^^ Mpc at z — 3). 

We do not know the redshift dependence of hrT' so we 
have to consider a range of possibilities and determine 
which leads to the most serious contamination. This is 
done for power laws bxT' cx (1 -I- z)^ in the last part of 
Table IIIII (note that fiat spectrum sources with constant 
comoving density would have C, — 2). Note that the 2b- 
and 3-halo terms are shown to be negligible. For C = 
we find upper limits to bTT'{x) to be (0 ± 5)/i nK Mpc^^ 
for the LRGs, (— 6± ll)ft,nKMpc~^ for the quasars, and 
(-1-1 ± 2)/inKMpc"^ for NVSS; these are upper limits 
because they assume the entire correlation comes from 
the 2-halo term. The combined 95% confidence upper 
limit is 4ft,nKMpc~^. This is a one-tailed upper limit 
{+1.64(t) since physially bxT' should be positive. 

We now apply Method II f Sect ion IVirXS)) to the l+2a 
halo term. For the method based on Ka — V and Q — W 
differences, and for fiat-spectrum sources the flux ra- 
tios {FKa - Fv)/Fv and {Fq - Fw)/Fv are 2.20 and 
1.63, respectively, so that AT^^^AT^^) = 3.57. The trans- 
fer function T/ varies from 5.0 x 10'^ to 3.0 x 10^ as I 
varies from to 400. The implied bias A A is given in 
Table IIIII We have shown this for both cases where 
the Galactic foreground model is subtracted ( "Gs" ) and 
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where it is removed by high-pass filtering ("hpf"). We 
have also shown this for different combinations of fre- 
quency bands. In cases where the same band appears 
twice, e.g. (Ka— V)(Q— V), we have computed only cross- 
correlations involving different DAs and averaged the 
results, e.g. "(Ka— V)(Q— V)" is actually the average 
of (Ka-V1)(Q-V2) and (Ka-V2)(Q-Vl). Just as for 
Method I we compute error bars using simulated LSS 
maps since it is harder to simulate the product-difference 
maps. 



these low-redshift sources must be constrained in order 
to have a reliable result. 

Note that many dusty sources at z > 0.02 will appear 
in the IRAS maps (in the sense that the IRAS maps 
contain their statistical fluctuations, even if the sources 
themselves are confused or buried under Galactic emis- 
sion). This means that these sources are double-counted 
in our argument, once in the statistical investigation of 
Cf and once in the foreground tests using IRAS maps. 
This does not represent a problem since it only serves to 
make our bounds more conservative. 



C. Estimates of contamination: infrared sources 

Infrared sources present a different challenge than ra- 
dio sources. On the one hand, our direct constraints on 
^max are quite weak due to the lack of coverage at fre- 
quencies above W band: there are no all-sky maps at 
frequencies between the WMAP W band (94 GHz) and 
COBE/DIRBE channel 10 (1250 GHz). However most 
infrared sources are quite weak (there are no sources in 
the WMAP catalog with spectral indices consistent with 
thermal dust emission) so taking i^max from the WMAP 
detection threshold would be very conservative (in fact, it 
would not yield a useful constraint). Therefore Method I 
by itself is not very helpful. On the other hand, the same 
low flux of the point sources makes Method II useless as 
any frequency-difference maps are dominated by Galac- 
tic emission on large scales and radio sources on small 
scales. 

In the face of these difficulties, we have resorted to a 
combination of two different approaches. The key idea 
behind our approach is that infrared sources at cosmo- 
logical distances have fluxes far below the WMAP de- 
tectability threshold. (Note that this is very different 
from the behavior of radio AGNs.) Therefore we con- 
sider separately the distant sources (z > 0.02) and nearby 
sources (z < 0.02). The distant sources are constrained 
using the arguments of Section IVII Al including Method 
I for the l+2a-halo terms. The redshift cut allows us to 
use a small value of i^max (essentially the flux of an object 
a few times L^, at z = 0.02). The nearby infrared sources 
can be handled by a very different argument: they appear 
in the IRAS 100 /xm maps, so we can use these maps to 
constrain them, i.e. we can feed the IRAS maps through 
our lensing pipeline with an appropriate scaling factor. 
(Distant infrared sources at z ~ 1, which are a major 
source of concern for us, may appear in WMAP but not 
IRAS because the fc-correction for a thermal dust spec- 
trum acts to brighten these sources in V and W bands, 
but make them fainter in the IRAS bands.) Being at 
different redshifts, the two groups of sources cannot have 
any physical correlation with each other. The sources at 
z < 0.02 should have very little overlap with our large- 
scale structure samples, especially the LRGs, but there 
are a few quasars with low spectroscopic redshifts, and 
NVSS must have a significant tail at z < 0.1 because of 
the nonzero NVSSx2MASS cross-correlations. Therefore 



Sources at z > 0.02 



We now apply the methods of Section IVII Al to the 
infrared sources. The infrared sources are subdominant 
contributions to the power spectrum in all of the WMAP 
bands, so we cannot easily implement Method II. On the 
other hand, they are much fainter than the radio sources, 
so Method I (based on -Fmax) is extremely useful. We will 
reference our constraints on their contamination to the 
W band since infrared sources are brightest there. For 
Fi, cx i/^-^ spectra and the relative weighting of bands in 
Eq. p7)) . the ratio of contamination is 



ACf (TT) 
ACf^WW) 



0.51; 



(66) 



for cx the coefficient is 0.44. 

Just as we did for the radio sources, we may construct 
a linear combination of differences Cf"^(Q) — Cf^(V^) and 
Cf^{V)-Cf^iW) to estimate Cf^ (W) without any con- 
tamination by the ISW effect. The difference between V 
and W bands is the most sensitive to infrared sources, but 
it contains an opposite contribution from radio sources, 
which must be canceled by involving the Q band. The 
appropriate combination is 

Cf(H^;irps) ^ a3[Cf (Q) ~ (!/)] 

+a4Cf^{V)-C^^{W)], (67) 

where 03 = -1-1.086 and 04 — —2.295. This combi- 
nation has no sensitivity to flat-spectrum sources and 
has the correct normalization for sources with spectrum 
Fi, cx i'^'^. If the spectrum of the infrared source is 
steeper than v^-^ (e.g. v"^) then this equation overes- 
timates Cf "^(VF; irps). If there is a positive contribution 
to Cj""^ from free-free or steep-spectrum sources {F^, cx 
with a < 0) then it adds a positive contribution to the 
estimate C'j''"^(M^; irps). Also it is easily verified that for 
a tSZ spectrum with blackbody temperatures in the ra- 
tio Tq : Tv : Tw = 0.957 : 0.906 : 0.783 and negative 
amplitude {Cf'^ < as appropriate since tSZ produces a 
decrement in high-density regions), a positive bias is pro- 
duced in Cf^ {W; irps). Therefore Eq. ([67]) represents an 
upper limit on the contribution to Cf'^ (W) from infrared 
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TABLE III: Estimates of the radio point source contamination \AA\ to the galaxy-convergence correlation amplitude A. The 
estimates are referenced to the VV lensing maps; for flat-spectrum sources the contamination to TT is reduced by a factor 
of 0.56. We consider separately the contribution from the l-|-2a, 2b, and 3-halo terms. The first part of the table consists of 
upper limits; the subscript "max" is used to emphasize that these are estimates of an upper limit. In particular |A^|max can 
be used to assess the magnitude of systematic errors but should not be used to "correct" the data! The bottom part of the 
table contains estimates of the contamination - these are "central" values and their errors (where shown) are la statistical 
errors, not including any uncertainty in the frequency dependence of the foreground. The notations "Gs" and "hpf" are used 
to indicate that the Galactic foreground is subtracted (Gs) from the difference maps D^^''^\ or that it is removed by high-pass 
filtering (hpf); see Section rVlTXSl 



Term 


Method/assumption 


lAAUax(LRG) 


|Ayl„,ax|(QSO) 


|AA|„,ax(NVSS) 


2b halo 


( = -1-4, fiducial cosmology 


8.8 X 10"* 


2.8 X 10~^ 


4.0 X 10-^ 


2b halo 


( ~ 0, fiducial cosmology 


5.6 X 10"^ 


6.2 X 10"^ 


7.0 X 10"^ 


2b halo 


( — —4, fiducial cosmology 


1.8 X 10"^ 


1.6 X 10"^ 


1.9 X 10"^ 


3 halo 


^ = 0, fiducial cosmology 


1.2 X 10"^ 


1.4 X 10"* 


4.3 X 10"* 


3 halo 


( — +4, fiducial cosmology 


1.7 X lO"'' 


6.5 X 10-* 


4.2 X 10"* 


3 halo 


( — —4, fiducial cosmology 


1.7 X 10"^ 


1.4 X 10"* 


1.1 X 10"^ 


Term 


Method/assumption 


Ayl(LRG) 


AA(QSO) 


Ayl(NVSS) 


l-|-2a halo 


Method II, (Ka-V)(Q-W), hpf 


+0.060 ± 0.045 


+0.039 ± 0.040 


+0.041 + 0.037 


l-|-2a halo 


Method II, (Ka-V)(Q-W), Gs 


+0.001 ± 0.076 


-0.007 + 0.062 


-0.023 ± 0.034 


1-f 2a halo 


Method II, (Ka-V)(Q-V), hpf 


+0.005 


-0.013 


+0.046 


l-|-2a halo 


Method II, (K-Q)(Ka-V), hpf 


+0.032 


+0.013 


+0.021 


l+2a halo 


Method II, (Ka-V)(Q-V), Gs 


-0.090 


-0.009 


-0.026 



We use Method I here for the l+2a-halo terms, so 
we need an estimate of i^max- For the radio sources we 
used the detection threshold in the WMAP point source 
catalog, however this hmit is too weak to be useful for 
constraining infrared sources because WMAP does not 
have any bands at higher frequency than W and conse- 
quently the sensitivity to infrared sources is quite weak. 
(Indeed, the WMAP point source catalog [ll contains 
no objects consistent with being infrared rather than ra- 
dio sources.) Therefore we take a different approach. 
The brightest sources (in Jy) may be either local ob- 
jects (bright because they are nearby) or high-redshift 
submillimeter galaxies. For the local objects we use the 
local SCUBA luminosity function at 353 GHz (the lowest 
measured frequency) from Dunne et al. [s^, which im- 
plies that within our survey solid angle and at z > 0.02 
there should be on average 0.05 objects with F > 0.5 Jy. 
Thus we take Fmax = 0.5 Jy at 353 GHz, which is ex- 
ceeded with 5% probability. High-redshift submillime- 
ter galaxies have typical fluxes of up to tens of mJy at 
353 GHz [H, [qO] and hence do not increase our above 
estimate of i^max- In order to be useful, this estimate of 
the maximum luminosity of the infrared sources must be 
extrapolated down to W band. Unfortunately there is 
very little data to suggest the correct form of the SED 
in this frequency range. This depends on the value of 
the emissivity exponent /? for the infrared sources (recall 
that Fi, (X v^^^ in the Rayleigh- Jeans limit). While low 
values of /3 ^ 1.3 are derived by Dunne et al. [SSJ from 
the SED fits for local galaxies (recall that F^ oc v'^'^^ 
in the Rayleigh- Jeans limit), Dunne et al. also note that 
these values could be biased low by a superposition of sev- 



eral dust temperatures, and that fits to our own Galactic 
emission give larger (3. We have assumed /3 = 1.3 to 
be conservative, which gives Fmax = 0.013 Jy in W band 
(this number decreases if (3 is increased). In blackbody 
temperature units this is i^max — 0.06 nK sr. 

The constraints we obtain are shown in Table HVl The 
l+2a and 2b-halo terms are based on the 95% confi- 
dence upper limit to — with the specified value of 
-Fmax (since AA is negative). 

For the 3-halo terms we need bxT' . Unfortunately in 
this case the galaxy-temperature spectrum is not good 
enough to constrain bxT' because for infrared sources 
bxT' is probably a strongly increasing function of red- 
shift. This can be seen from Eq. if the comoving 
density of infrared emitters is constant, then due to the 
/c-correction we have T" cx (1 + z)^+" where the spectral 
index a is typically '--^ 3.5. Also the bias is usually higher 
for high-redshift galaxies where the fluctuation amplitude 
is less and hence galaxies can only form in rare, highly bi- 
ased peaks of the density field. Finally the star formation 
rate and hence (probably) the infrared luminosity density 
declines at z < 1 (e.g. f91]). Given the very large un- 
certainties in the redshift distribution of the far-IR back- 
ground, we have tried several models for T". In Model A, 
we assume that at 1 < z < 5 the comoving infrared emis- 
sivity is constant and a — 3.5 so that f ' cx (1 + z)**-^; 
at z < 1 it is assumed that the emissivity declines as 
oc (1+z)'^ so that T" c>c (1+z)^'^. We normalize the model 
to produce 33 /iK blackbody temperature at W band over 
the redshift range z < 5, consistent with extrapolation of 
the far-IR background observed by COBE/FIRAS [9^. 
This normalization may be conservative: there must be 
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some contribution to the far-IR background from z > 5, 
particularly at low frequencies where the fc-correction fa- 
vors high redshifts. The other models include several 
variations on this. Model B does not include the decline 
at low redshift, and takes T" c>c (1 + z)^-^ all the way 
down to z = 0; it is designed to maximize the contam- 
ination for the LRGs. Models C and D are similar to 
Model A except that the exponent of 1 -I- z is increased 
or decreased by ±2 respectively. For the bias, we have as- 
sumed a constant clustering strength with as{z, gal) = 2, 
i.e. b = 2/(TgG{z) where G is the growth factor. This 
is somewhat larger than what is observed for the LRGs 
and quasars, and is probably typical for the most mas- 
sive objects. The actual emissivity-weighted mean bias 
is probably somewhat lower. The resulting estimates of 
the 3-halo infrared source contamination are shown in 
Table HVl While the calculation is obviously very crude, 
it is clear from the table that the infrared source 3-halo 
term is not a significant source of contamination. 



2. Sources at z < 0.02/ constraints from IRAS 

The nearby infrared sources are constrained u sing the 
100 /im maps from the IRAS Sky Survey Atlas [Mlsl' 
which cover ~98% of the sky with ~6 arcmin resolution. 
Our strategy is to feed a rescaled version of the IRAS 
maps through our lensing pipeline, and use the result 
as an estimate of the contribution to Cf" from nearby 
dusty galaxies. (More accurately, we do the rescaling 
assuming a low dust temperature to be conservative, and 
treat this as an upper bound on the contamination.) We 
note that the IRAS maps are dominated by emission from 
our own galaxy: this should not represent a problem, as 
the output of the lensing pipeline should then contain 
the contamination to v from both Galactic and nearby- 
extragalactic dust. 

In order to implement the above program, we need 
to select an SED with which to extrapolate from 
100 /im to V/W bands, in particular we need the ratios 
Iu{V)/I^ilOO^im) and I^{W)/I^imOnm). If we assume 
that the nearby galaxies have SEDs similar to the high- 
latitude regions of the Milky Way, and use the "cold re- 
gion" spectrum from Fig. 1 Finkbeiner et al. [s^, we esti- 
mate I^{W)/I^{100fim) = 0.0012. The regions of higher 
dust temperature would have a smaller ratio (bluer spec- 
trum) and hence this is conservative. We mask Saturn 
from the IRAS maps since planet-contaminated data is 
not used in WMAP mapmaking. 

The resulting contamination estimates in W-band are 
AA = +0.002 for the LRGs, AA = +7 x 10"^ for the 
quasars, and = -f 0.009 for NVSS. This is negligible 
compared to our error bars. 



D. Thermal SZ effect 



The thermal SZ effect is another possible contaminant 
of the lensing signal. This section evaluates contamina- 
tion from the galaxy-tSZ-tSZ bispectrum. The method 
of setting constraints for the 2b- and 3-halo terms is very 
similar to the approach for point sources, except that 
we estimate bxT' from theory rather than from cross- 
correlation of the galaxies with the temperature field. 
This is because on the one hand the weak frequency 
dependence of tSZ (within the WMAP bands) makes 
the observational measurement of the galaxy-tSZ cross- 
correlation very noisy; and on the other hand we have 
at least a rough idea of how and why the tSZ decre- 
ment depends on halo properties (as opposed to point 
sources, whose halo occupation properties can only be 
determined empirically). The 1- and 2a-halo terms can 
also be constrained theoretically from the mass function 
and mass-bias relation (for the 1-halo term we also need 
Nniax)- They diverge at low redshift and so we have esti- 
mated the contribution from nearby haloes by using their 
actual positions from X-ray data. The bounds in this sec- 
tion will be written in terms of the Rayleigh- Jeans (RJ) 
limit of the tSZ effect; since WMAP is not quite in the 
RJ limit, the galaxy-tSZ-tSZ bispectrum is suppressed 
by a factor of 0.82 (for the VV frequency combination), 
0.71 (VW), 0.61 (WW), or 0.72 (for TT, the frequency- 
averaged lensing map). 

The nearby haloes are obtained from the Ebeling et al. 
XBAC catalog 95] of optical clusters detected in the 
ROSAT 0.1-2.4 keV X-ray band. We have used the X-ray 
temperature to SZ luminosity conversion normalized ac- 
cording to Hinshaw et al. [l| [i.e. we used their Eq. (29) 
multiplied by either the central value of the normaliza- 
tion —0.32, or the one-sided 95% confidence lower limit 
—0.55]. Some of the redshifts in the XBAC catalog are 
proprietary but there are now publicly available redshifts 
for these clusters [ol, [93, [11] , which we use. With this 
data it is possible to construct the actual realization of 
the tSZ sky contributed by the nearby X-ray luminous 
clusters, and to feed this through the lensing pipeline 
to obtain the contamination AA, if one knows what as- 
sumption to make for the tSZ radial profile of the clus- 
ters. We found earlier that Gaussian or NEW profiles 
always produce smaller lensing contamination than delta 
functions, so it is conservative to treat the cluster as a 
delta function. We do this for all clusters except the 
brightest (Coma); following Hinshaw et al. [l| we have 
modeled Coma using a (3 profile [99] with central tem- 
perature decrement at 94 GHz of either —0.4 2 mK (best 
value) or -0.54 mK (95%CL, one-sided) poj . 

The 2b- and 3-halo terms for the tSZ effect depend on 
the product bxT' for the tSZ effect. In the RJ limit, and 
assuming that the intracluster gas contains the cosmic 
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TABLE IV: Estimates of the infrared point source contamination to the galaxy-convergence correlation amplitude A. The 
estimates are referenced to the WW lensing maps; for cx sources the contamination to TT is reduced by a factor of 0.51. 
We consider separately the contribution from the l-|-2a, 2b, and 3-halo terms. The nearby sources are not accounted for in our 
analysis of the l-|-2a halo terms, so we have used the IRAS maps for these; those estimates are shown in the last line of the 
table, and are estimates of A^, not upper limits. 



Term 


Method/ assumption 


|A^|„,ax(LRG) 


lAAnaxKQSO) 


|A4|,^ax(NVSS) 


l+2a halo, z > 0.02 


Method I 


0.0089 


0.0087 


0.0003 


2b halo 




0.0024 


0.0062 


0.0040 


3 halo 


Model A, fid. cosmology 


6.6 X 10"^ 


6.0 X 10"^ 


1.2 X 10"* 




Model B, fid. cosmology 


3.3 X 10"® 


6.0 X 10"^ 


1.2 X 10"* 




Model C, fid. cosmology 


6.6 X 10"^ 


7.3 X 10"^ 


9.1 X 10"^ 




Model D, fid. cosmology 


4.6 X lO"'^ 


4.4 X 10"* 


3.1 X 10"* 




Model A, high ag 


8.7 X lO"'^ 


6.6 X 10"^ 


1.3 X 10"* 




Model A, high f7„ 


1.6 X lO"*' 


1.0 X 10~* 


1.8 X 10~* 


Term 


Method/ assumption 


Ayl(LRG) 


Ayl(QSO) 


AA(NVSS) 


sources at z < 0.02 


Extrapolated from IRAS lOO/^m map 


+2 X 10~^ 


+7 X lO^'""' 


+9 X lO^-^ 



baryon fraction, this is given by 



fb 



irieC 



(68) 

where r' is the mean optical depth per unit comoving dis- 
tance, h{M) is the bias for haloes of mass M, fee is Boltz- 
mann's constant, me is the electron mass, Te{M) is the 
density-weighted mean electron temperature in haloes of 
mass M , and (piM) dM is the fraction of the mass in 
haloes of mass between M and M + dM. Here ficu/fb 
is the ratio of the intracluster medium baryon fraction to 
the cosmic value. This is < 1 in real clusters, since some 
of the gas has turned into stars or is present in atomic or 
molecular phases that do not contribute to tSZ; however 
we will set ficu/fb = 1 here to be conservative. Note 
that / (j){M)dM = 1. 

In evaluating this expression we use the Jenkins et al. 
[lOll | mass function (fit to A^-body si mula tions) and the 
Sheth & Tormen mass-bias relation [TO^ . The density- 
weighted gas temperature cannot be obtained from A^- 
body results; we have used the analytic model of Reid & 



Spergel [Tosj for the SZ luminosity, which after applica- 
tion of their Eq. (2) yields 



n = 5.9 



10i5/i-iM( 



o 



2/3 



M 



keV. 



(69) 



Note that here we use the virial mass My based on 
the spherical collapse model for consistency with Reid 
& Spergel 103], whereas the Jenkins et al. [Toii mass 
function uses M defined based on a spherical overden- 
sity of 180p. We have used the NFW profile and mass- 
concentration relation [l04| to convert M to My. The 
fraction My/M is nec essa ry since the total SZ luminos- 
ity predicted by Ref. jl03l | is actually the product MTe. 
The corresponding bounds from Section FVII Al are shown 
in Table El 

We now turn to the evaluation of the 1- and 2a-halo 
contamination, from Eq. (|36p . The 2a- halo term can ac- 
tually be calculated provided we know the background 
cosmology, mass function, bias, and mass-temperature 
relation: 



ACr(2ah) = J dFj2^n2D{N)n2D{F)QiN;F)^^jj^ J dx J dF F\n{x)nMF)b{F)PUk) 



rpsjl ) f dx 

Ri 



I % [ dM^m)KM)fix)PUk) (2f'T,uB^'-^I^y , 
J r'^ J Po V Jb meC^ J 



(70) 



where the first equality involves use of the Limber ap- 
proximation and the fact that the (2-halo) galaxy-halo 
cross spectrum is the product of biases times the linear 
power spectrum, and the second equality involves conver- 



sion of the integral from tSZ flux F to halo mass M. The 
denominator po is the mean matter density of the uni- 
verse today. The argument for the 1-halo term is similar 
except that here we actually need the number of galaxies 
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per halo, and there is no factor of Ci{N; F): 



ACT(lh) = 



rps{l) f dx 



Ri 



M 



-1 / dM—(j){M) 
po 



N{M) 



^(.y-lrr. /iCM kjiT^{M) \^ 

X 2t JcMB— J 5 — ■ (71j 

V h meC-^ J 

We may evaluate Eqs. (jTO]) and (ffTj) in a worst-case sce- 
nario by assuming ficu/fb = 1 and N{M) = A^,„ax- Not 
surprisingly, these equations are dominated by the low- 
redshift portion of the integral (z < 0.2), especially for 
NVSS where there is no photo-z cut (such cuts are im- 
possible for radio continuum-selected objects). Unfor- 
tunately these equations lead to rather weak constraints, 
especially on the 1-halo term. Therefore for the results in 
the table we multiply the integrand by 1 — PxbaC; where 
^XBAC is the probability that a cluster with a given mass 
would pass the XBAC flux threshold. This gives us the 
1- and 2a-halo contributions from clusters whose X-ray 
flux is too low to have been included in XBAC; conse- 
quently our estimates of contamination should be added 
to those obtained from XBAC. The probability Pxbac 
was then determined using the XBAC flux li mit a nd the 
lognormal P{L\M2oo) distribution from Ref. [l05j |. 



E. Summary of extragalactic foregrounds 

The analysis of extragalactic foregrounds presented 
here has been quite long, and it is worth summarizing 
our major findings. These have been that the 1, 2a, and 
2b-halo terms for the radio sources and tSZ effect could 
contribute at the level of up to several percent. For the 
infrared sources, the contamination \AA\ (rescaled to the 
frequency-averaged map) is < 1% for each of the samples, 
even after adding 1, 2a, 2b, and 3-halo terms and the local 
contribution estimated from the IRAS maps. The radio 
and tSZ 3-halo terms are smaller (< 0.2%, after rescal- 
ing to the contamination in the cross-correlation with the 
frequency-averaged TT lensing map); and the kSZ and 
Rees-Sciama effects are completely negligible (< 0.1%). 
This is good news, since the estimates we have used for 
the 3-halo terms and the kSZ and RS effects are crude 
(based on perturbation theory in the quasilinear regime) 
and we would not trust the theory to give a correction 
for them. 

We now discuss the potentially significant foreground 
terms and the robustness of their contamination esti- 
mates: 

1. For the radio sources, the contamination may be 
dominated by the l-|-2a halo term (the uncertainty 
in this term is larger than the maximum effect from 
the 2b-halo term). These have been constrained 
by constructing a product of frequency-difference 
maps D (which contains a feature around each 
point source whose amplitude is proportional to the 
square of the flux) and cross-correlating this with 



the galaxies. There is no detection of this correla- 
tion Cf^: the values shown in the bottom section 
Table IIIII are consistent with zero. The 2a worst- 
case contamination is for the LRGs where \AA\ 
could be as large as 0.15. 

2. For the tSZ effect from nearby clusters, the con- 
tamination was assessed primarily using the XBAC 
catalog. The estimate for the former was normal- 
ized using the WMAP results Q and the result is 
that the contamination is at the \AA\ < 0.01 level 
in the frequency-averaged maps. 

3. The tSZ effect from distant or faint clusters below 
the XBAC flux limit was constrained using theo- 
retical arguments: for a given cosmology we can 
calculate the halo mass function and bias. The 
worst contamination is for the high-tJs cosmology 
(cTg — 0.92) so we focus on this case. The 1-halo 
term in Table fVl is the largest (up to 3%), how- 
ever this is an upper limit assuming that there are 
NmsLX galaxies in every massive cluster. This may 
not be far from the case for NVSS, however even 
for massive clusters the LRG count is typica lly a 
factor of several smaller than A^max = 36 jl06l | 
and the photo-z cuts should remove most LRGs 
from z < 0.2 clusters. Also most clusters will host 
much fewer than A^max — 10 quasars - although 
one should remember that occasionally non-quasar 
extragalactic objects get counted as quasars due 
to failures in the photometric pipeline (e.g. the 
pipeline shredded the galaxy NGC4395 and some 
of the Hii regions were classified as quasars). In 
contrast, the 2a-halo term is a best estimate: the 
contamination is ~ 1.4% for the LRGs and less for 
the quasars and NVSS since they live at higher red- 
shift where there are many fewer massive clusters. 
While it has some uncertainty due to the Lsz{M) 
relation, it does not dominate and so we believe the 
sum of all the terms (1, 2a, 2b, 3-halo, and XBAC) 
provides a reasonable upper bound on tSZ contam- 
ination. This is \AA\ < 0.057 for the NVSS-TT 
correlation, and less for the LRGs and quasars. 

4. The arguments in Appendix [Dl show that the cross- 
correlations of different foregrounds are not impor- 
tant since none of the foregrounds are important 
individually. 



VIII. DISCUSSION 

In this paper, we have measured the cross-correlation 
between three samples of galaxies (LRGs, quasars, and 
radio sources) and lensing of the CMB. We find evidence 
for such a correlation at the 2.5ct level, with an am- 
plitude of 1.06 ± 0.42 times the expected signal for the 
WMAP cosmology. All sources of systematic error and 
foreground contamination are believed to be negligible 
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TABLE V: Upper limits to the contamination to the galaxy-convergence spectrum from the tSZ effect. The 2b- and 3-halo 
terms have been evaluated for the three cosmologies considered. The contamination estimates are scaled to the Rayleight- Jeans 
limit; multiply by 0.72 to get contamination to the frequency- averaged lensing signal. The local XBAC cluster constraints are 
computed using several assumptions about the radial profile, and several assumptions about the fluxes: pessimistic ("pess", 
based on 95% confidence upper limits) and central values ("cent"). 



Term 


Method 




AA(LRG) 


AA(QSO) 


Ayl(NVSS) 


2a halo 


fid. COS. 




-0.0024 


-0.0001 


-0.0010 


(nonlocal) 


high as 




-0.0202 


-0.0008 


-0.0061 




high Qm 




-0.0030 


-0.0001 


-0.0013 


Local clusters 


XBAC-cent, /? (Coma) 


1, pointlike (all others) 


-0.0007 


+0.0007 


-0.0039 


Local clusters 


XBAC-pess, f3 (Coma) 


1, pointlike (all others) 


-0.0045 


+0.0032 


-0.0146 


Term 


Method 




|Ayll„,ax(LRG) 


lAAUax(QSO) 


lAA|„,ax(NVSS) 


1 halo 


fid. COS. 




8.f X fO~^ 


7.5 X 10"^ 


8.7 X 10"^ 


(nonlocal) 


high (jg 




4.f X fO"^ 


3.8 X 10"^ 


4.4 X 10"^ 




high Qm 




f .6 X fO"^ 


1.5 X 10~^ 


1.7 X 10"^ 


2b halo 


fid. COS. 




1.3 X 10"^ 


1.4 X 10~^ 


1.6 X 10"^ 




high as 




1.2 X 10"^ 


1.2 X 10"^ 


1.4 X 10"^ 




high 




1.1 X 10"^ 


1.0 X 10"^ 


1.2 X 10"^ 


3 halo 


fid. COS. 




3.1 X 10"^ 


2.9 X 10"^ 


1.0 X 10"* 




high as 




3.0 X 10"^ 


2.6 X 10"* 


1.0 X 10"^ 




high Qm 




8.7 X 10"* 


4.9 X 10"^ 


6.0 X 10"* 



compared to the statistical uncertainties. Our measure- 
ment is consistent with the earlier nondetection by Hirata 
et al. [23] using LRGs and the more recent 3.4ct result 
by Smith et al. [Ill using radio sources. 

The case for lensing of the CMB is strengthened by 
having analyses by two different groups (Smith et al. and 
us). While the results both draw on WMAP and NVSS 
data, and hence are not independent, it is worth noting 
the very significant differences between the analysis pro- 
cedures. In the lensing reconstruction procedures, Smith 
et al. used weighting of the WMAP data, whereas 
we force the reconstruction to use the same weight func- 
tion Wi at all frequencies and in regions with different 
noise levels in order to simplify the foreground analy- 
ses. They also use autocorrelations of the same CMB 
map and subtract off the noise-induced biases, whereas 
we use cross-correlations which inherently have no noise 
bias. Our NVSS sample was constructed with a flux cut 
at 2.5 mJy whereas Smith et al. use the full catalog (ex- 
cept for masked regions) . We also used cross-correlations 
with other galaxy samples to estimate the NVSS bias and 
redshift distribution (see Paper I) , instead of fitting the 
auto-power spectrum; such cross-correlations are usually 
more robust, especially for a sample such as NVSS that 
shows large amounts of instrumental power. While even 
large uncertainties in the redshift distribution do not al- 
ter the detection significance (in number of sigmas) , with- 
out a solution to the redshift distribution problem there 
is no hope for CMB lensing in cross-correlation to ever 
become a useful cosmological probe. 

We would also like to test for consistency between the 
two results. Smith et al. found a cross-correlation of 



1.15 ± 0.34 times the expected signal. However, because 
their redshift distribution model is peaked at higher z, 
and they have a different fiducial cosmology (higher cts), 
their predicted signal amplitude is larger than ours; fit- 
ting the Smith et al. cross-spectra Cf^ to our theoretical 
prediction gives an amplitude oi A — 1.32 ± 0.40. This 
compares well with our result of Anvss = l-H ± 0.52. 
[The results are still not quite comparable because of the 
flux cut difference (2.5 versus 2.0 mJy), although we sus- 
pect this would not cause a huge difference because only 
18% of the NVSS catalog is at < 2.5 mJy, and the bias 
and redshift distribution of radio sources is usually be- 
lieved to change very slowly with flux [sij . In particular 
to explain the entire difference in amplitude, the lens- 
ing cross-correlation signal A from the faint F < 2.5 mJy 
sources would have to be a factor of ~ 2 higher than 
for the bright sources, which we consider unlikely.] The 
results were also arrived at independently: no revisions 
have been made to our central value since we became 
aware of the Smith et al. result (although our error bars 
were not finalized at that time because we had not com- 
pleted analyzing our simulations). An obvious question 
in the comparison is why the Smith et al. error bar is 
smaller than ours by a factor of 1.3. Some of this is due to 
the 20% higher number density in their Poisson-limited 
NVSS map due to their lower flux cut, which accounts 
for a factor of \/l-2. The remaining difference - a fac- 
tor of 1.2 in standard deviation - must be attributable 
to all other differences in the analysis, such as their use 
of Q band and auto-correlations between WMAP maps, 
different handling of point sources, and (going the other 
direction) Smith et al. inclusion of the point source and 



26 



Galactic foreground systematics (which we find in this 
paper to be negUgible). 

Because of the different choices made in the analysis 
the final signal amplitudes end up being different and 
the Smith et al. central value is higher. This must be 
a statistical fluctuation if both methods are unbiased, 
unless the bdN/dz is dramatically larger for the fainter 
sources. These sorts of fluctuations can matter if one is 
trying to claim a detection, e.g. our measured ampli- 
tude would only give a 2.8ct signal even with their errors 
versus their 3.4tT. However, in this paper we take the 
view that what is ultimately relevant is the power of the 
method to disciminate between the cosmological models 
and the advantage of having a higher detected amplitude 
can quickly be turned around if there is an interesting 
class of models that predicts a larger lensing signal than 
the standard model, and which would therefore be more 
strongly ruled out if the measured amplitude is lower. It 
is for this reason that we put such an effort in determining 
redshift distribution weighted bias bdN/dz for the sam- 
ples, because without it no cosmological interpretation 
is possible, since variations in cosmology and in b dN/dz 
will be degenerate. 

Another difference in the two analyses is in the treat- 
ment of extragalactic foregrounds. We split the fore- 
grounds into 1, 2a, 2b, and 3-halo terms for each fore- 
ground component, and considered each of the many re- 
sulting terms separately. Smith et al. did several tests, 
the most important being (i) a comparison of Q, V, and 
W band signals, and (ii) a point source bispectrum test 
which in our terminology accounts for the l+2a halo 
terms. There are possible worries with these tests, e.g. 
in their implementation of (ii) the "template" bispectrum 
for the point sources has the same frequency dependence 
as the CMB. Real point sources will have a different fre- 
quency dependence, however the template should still 
have some sensitivity to them and this is demonstrated 
explicitly by the Smith et al. simulations. The method 
also does not consider 2b and 3-halo terms, and both 
tests can miss kSZ or RS contaminants; however as we 
have argued here these are negligible. 

We also included the SDSS LRGs and quasars. These 
together (without NVSS) give an amplitude of 0.99±0.56 
times the expected signal; this is a l.Scr result and is 
consistent with the fiducial cosmology. Thus in all cases 
the measured lensing amplitude is consistent with the 
WMAP cosmology, which can be counted as a success. 
However it is not yet competitive with other cosmological 
probes. We have constructed a likelihood function (see 
Appendix |E| and included it in the Markov chains of Pa- 
per I as a proof of principle, but it does not add much to 
the CMB-I-ISW constraints. The real challenge for lens- 
ing of the CMB is to move beyond the "first detection" 
stage to providing interesting cosmological constraints. 
Several improvements must be made in order to bring 
lensing of the CMB to the leading edge of cosmology: 

1. The most obvious issue is that CMB maps with 
higher resolution and lower noise than WMAP are 



required (one cannot improve on the sky cover- 
age). Th is wi ll get better with the upcoming Planck 
satellite (ll4j and the g round-based Atacama Cos- 
mology Telescope jllSl ] and South Pole Telescope 

[Hi. 

2. Extragalactic foregrounds, while only a minor issue 
here, get worse as one uses the higher multipoles 
in the CMB. One possibility is the use of multi- 
frequency information: at the higher frequencies 
probed by the future experiments, infrared point 
sources and tSZ have a distinctive spectral signa- 
ture. Another possibility is to use spatial infor- 
mation: a modified reconstruction procedure could 
be used to suppress the point source responsivity 
function rps(Z), which controls the 1- and 2a-halo 
contributions from point sources. (It is however im- 
possible to do a reconstruction that is insensitive 
to all possible forms of the 3-halo term.) Yet an- 
other possibility, which is probably the best when 
the data become available, is to use polarization 
[m, [ill. The point sources, tSZ, and kSZ effects 
are either observed or theoretically expected to be 
m uch weake r in polarization than in temperature 
[3,[io3,[ioi, and the RS effect is absent. Also for 
randomly oriented point sources the 2b- and 3-halo 
contamination is exactly zero because two different 
sources can have no QQ or UU correlations, thus 
(aside from correlations of position angles of dis- 
tinct sources) eliminating rps{l) or using the multi- 
template bispectrum estimator of Smith et al. |28| 
would remove the point sources entirely. 

3. The third issue for lensing in cross-correlation is the 
theoretical uncertainty in the galaxy-convergence 
cross-spectrum. In this paper we have used linear 
biasing, Pgs{k) — bgPgg^{k), and thrown out small- 
scale information (fc > 0.1/iMpc~^) where linear 
biasing is not valid. This must be improved in the 
future: especially for low-redshift mass tracers a 
large fraction of the information is in the nonlinear 
regime, and if one aims for constraints at the several 
percent level, simple ideas like linear biasing may 
not be adequate even at fc = 0.1/iMpc~^. There is 
also uncertainty in the redshift distribution of the 
galaxies. One obvious method is to do lensing of the 
CMB in a utocorrelation (with the four-point func- 
tion |l09| or using the iterative approaches 73]), 
although this will be very demanding on control of 
instrumental systematics. If one chooses to go the 
cross-correlation route, then when sufficiently high 
signal-to-noise ratio is available it will be desirable 
to use halo modeling, and to better calibrate the 
photometric redshifts for the large scale structure 
tracers with spectroscopic surveys. The latter will 
of course also be valuable to any cosmic shear pro- 
gram using the galaxies as lensing sources. 

If these challenges can be met, weak lensing of the 
CMB will make the transition from being a simple con- 
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sistency check of the cosmological model to a routinely 
used cosmological probe. 
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APPENDIX A: CONTAMINATION FROM 
EXTRAGALACTIC MICROWAVE SOURCES 

The purpose of this appendix is to examine how ex- 
tragalactic sources (radio or infrared point sources, or 
Sunyaev-Zel'dovich effect) affect the measurement of the 
galaxy-convergence cross-power Cf^ . 

From Eqs. PH)) and we can see that after averag- 
ing over the WMAP instrument noise we have contami- 
nation 



(-1)'"(2/ + 1) J2 -^"'^ 



I I' I" 



Tl'm'Tl"m" , (Al) 

where we have introduced the coupling coefficient 



The three-mode correlation function of the galaxies 
and the temperature is defined via the relation 



{gimTl'm'Tl" 



I V I" 



(A3) 



Combining with Eq. (jAl[) . and using the reality condition 
= we find 
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(A4) 



where in the last line we have used the normalization of 
the 3j symbol to collapse the sum over azimuthal quan- 
tum numbers. Therefore any bias produced by point 
sources in the bispectrum Bf^j^, produces a correspond- 
ing bias in the galaxy-convergence cross-correlation. This 
is expressed in Eq. 
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APPENDIX B: DERIVATION OF TERMS IN 
POINT SOURCE BISPECTRUM 

This appendix derives the various multi-halo terms in 
the galaxy-point source-point source bispectrum, used in 
Section IVII Al In each case, the starting point for the 
derivation is Eq. (jA3[) . 



1. 1-halo term 

The 1-halo term is the simplest to derive. If a halo 
containing N galaxies with flux F is found at position 
n, then the contribution to gim is {N/ng)Yi^{n) and the 
contribution to Tim is FF;J^(n). The 1-halo contribution 
to the 3-mode correlation function is then obtained by 
integrating (or summing) over F, N, and n: 



{gimTl'm'Tl 



l"m" /Ih 



= [dFj2^F\,,, 



N 



I V I'' 



m m m 



(Bl) 



(The usual expression in terms of Legendre polynomials 
of cos 9 is equivalent to this by the spherical harmonic 
addition theorem.) Then we have 

{gimTl'm'Tl"m")2£ih — ^ / dF^— 
LM •' N ^ 

x[4tt6lo + Cl{N;F)] 
X J d'hYLM{n)YCm{n) 

X J d?h'YlM{i^') 

xy;„,(fi')r;,,„„(n'). (B4) 

The n integral is collapsed using spherical harmonic or- 
thonormality to SuSMnn f^nd the n' integral is a 3- 
harmonic integral. Therefore for Z 7^ we have 

/N 
dFy^—F^Ci{N;F) 



xmi'i" ' . (B5) 

\ m m m 



The evenness oi I + I' + I" for nonzero rju'i" and reality 
of the 3j symbol are used here, just as in Eq. (jBip . This 
equation proves Eq. 



where 77;;'/" is the combination of coefficients in Eq. ([32]). 
In the second equality we have used the three spherical 
harmonic integral formula, and removed the conjugation 
symbol because the 3j symbol is real. Comparison to 
Eq. jXS]) proves Eq. 



2. 2a- halo term 

The 2a-halo term involves two haloes, one with N 
galaxies, and one with flux F. These two types of 
halo have an angular power spectrum Ci{N;F) and 
correspondingly have an angular correlation function 
w{9\N;F). The probability density to find these two 
haloes at locations n (for the halo containing the galaxy) 
and n' (for the halo containing the point source) is then 
modulated by the factor 1 -I- 'w{9): 



N 



x^F'[l + w{0)]Y,Uh) 



xr;„,(n')l^r,„„(n') (B2) 



where 9 is the angle between n and n'. This factor can 
be re-written as 

l+w{9)^Y.^4TTSLo+CL{N;F)]YLM{n)YlM{h'). (B3) 

LM 



3. 2b-halo term 

Here we consider the gTT bispectrum introduced by 
correlations of a halo at n with N galaxies emitting flux 
F, and a halo at n' emitting flux F'. The angular power 
spectrum of these two haloes is Ci{N, F; F'). For this 
case, we have in analogy to Eq. (jB2p : 

{gimTi'm'Ti>,m")2bh = J J dF' J ^^'^ / ^^'n' 

N 

x^FF'[l + w{9)]YCmin) 

Hg 

xr;„,(n)y;,„„(n') 

+{l'm' ^ l"m"). (B6) 

The symmetrization in the last lie occurs because the 
halo at n could contribute to Ti'm' and that at n' could 
contribute to Tinm", or vice versa. Manipulations similar 
to those for the 2a- halo term then lead to Eq. ([46]) . 



APPENDIX C: SMALL EXTRAGALACTIC 
FOREGROUNDS 

1. Kinetic SZ effect 

The galaxy-kSZ-kSZ bispectrum is another potential 
contaminant of the CMB lensing signal, and is especially 
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worrying for future experiments because the kSZ signa l 
has the same frequency dependence as CMB lensing jllOj l . 
This section constructs a rough upper Umit to the kSZ 
contamination of the lensing signal, and finds that this 
contamination is negligible in WMAP. 

The kSZ temperature anisotropy is a line-of-sight in- 
tegral, 



kSZ 



Tc 



MB 



dxf'{l + 5bH,\\, (CI) 



where St, and ui, are the baryon density perturbations and 
velocity || denotes the line-of-sight component of velocity 
(positive toward the observer), and f' is the mean Thom- 
son optical depth per comoving radial distance. The ve- 
locity field will be irrotational (except for contributions 
on very small scales due to nonadiabatic effects) so the 
integral of velocity cancels along the line of sight. We 
will also make the approximation that Sb is equal to the 
matter density S, which should again be valid except on 
very small scales. (For comparison, the dominant con- 
tribution to the lensing signal comes from I ~ 400, or 
k = 0.4/1 Mpc~^ at the 25th percentile redshift of the 
LRGs.) In this case, the kSZ temperature increment 
becomes J dxf'Su\^ and the galaxy-kSZ-kSZ bispec- 
trum is given by 



<"(kSZ) = 



VU'l"T^MB f dX-; 



^r-n(x) 



xB{k,k',k"\g;SufSu\\). (C2) 



Here n(x) is the comoving distance distribution of the 
galaxies. 

In order to proceed we make a further assumption that 
the velocity field is perfectly coherent, i.e. that all of the 
power in U|| is at very small wavenumber and we can ob- 
tain the statistics of the kSZ signal by multiplying the 
statistics of the density field by the appropriate power 
of some large-scale velocity W|| and then averaging (this 
gives for statistics involving two copies of the kSZ 

signal; the factor of 1 /3 comes from the fact that the kSZ 
effect depends on only 1 of the 3 components of velocity) . 
This idealization is commonly made in predictions of the 
small-scale kSZ power spectrum and should be valid for 
predictions on scales smaller than the peak of the velocity 
power, i.e the maximum of kP{k) at fc ~ 0.05/i Mpc"^. 
Since larger scales are considered here, especially for the 
high-redshift part of the integral relevant to the quasars 
and NVSS, this assumption is probably not valid - in- 
stead the velocity field will become incoherent and the 
radial velocities of two widely separated regions will have 
only a small correlation. This is not a problem for us 
since the loss of coherence suppresses the kSZ signal, i.e. 
we maximize the possible kSZ signal by ignoring it. The 
reader should thus keep in mind that the galaxy-kSZ-kSZ 
bispectrum derived here may be a large overestimate. We 



have thus converted Eq. (|C2p into 



<r(kSZ) = r^u'vTl^B^ J ^r'^n(x) 



^ CMB -3^2 



xB{k,k',k"\g;S;S). 



(C3) 



The simplest way to calculate the bispectrum is to use 
second-order perturbation theory assuming linear bias for 
the galaxies, i.e. 



mi' I 



^CMB 



dx. 



;'2 



3c2 

x/(x)i?2p(fc,fc',fc"), 



(C4) 



where /(x) is the product of the bias and comoving dis- 
tance distribution. The calculation of the implied con- 
tamination to the lensing signal can be carried out just 
as for the point source 3-halo term (Eq. [63]) : one only 
makes the replacement 



brT' 



Wrms _/ rr 
-T 1< 



CMB- 



(C5) 



Just as in Eq. we take the absolute value of each 

term in the sum over I' and I" to avoid accidental can- 
cellation of terms due to the configuration dependence 
of the bispectrum (which for this extremely crude cal- 
culation we do not trust), and use linear theory to ob- 
tain Urms- For the fiducial cosmology, linear theory gives 
Wrms = 430kms~^ today. The resulting maximum con- 
tamination to the lensing signal is found to be |A^| < 
lO^'^ for each of the three cosmologies and each sample. 
The calculation presented here is very crude, nevertheless 
it suffices to establish that the kinetic Sunyaev-Zel'dovich 
effect has no significant impact on our results. 



2. ISW/Rees-Sciama effect 

Another foreground signal is the Rees-Sciama (RS) ef- 
fect, i.e. the sourcing of temperature fluctuations by 
changing gravitational potentials. This is physically the 
same phenomenon as the ISW effect, although the label 
"ISW" is usually used to refer to the linear-regime ef- 
fect; for brevity we will use the label RS here to refer to 
the effect on all scales. It arises at low redshift and thus 
can produce a galaxy-RS-RS bispectrum. This section 
considers an order-of-magnitude estimate of the galaxy- 
RS-RS bispectrum and the consequent contamination to 
the lensing signal, and shows that it is negligible. 

The RS effect can be written as 



T. 



CMB 



(C6) 



where 4' is the dimensionless gravitational potential and 
the dot denotes derivative with respect to conformal 
time. It is related to the density field via the Poisson 
equation. 



*(k) 



2ac2fc2 



5(k). 



(C7) 
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For linearly biased galaxies, the galaxy-RS-RS bispec- 
trum can then be written as 



^f{x)Bik,k\k"\6-^;i'). 

(C8) 

The problem is then to estimate the 3-dimensional bis- 
pectrum fc', fc"|i5; Vt; ^f). All bispectra are zero in 
linear perturbation theory with Gaussian initial condi- 
tions, so we have used the tree-level second-order bis- 
pectrum. The contamination to the galaxy-convergence 
cross-spectrum can then be estimated using a method 
analogous to Eq. 



(C9) 



As before, the absolute values are used to prevent acci- 
dental cancellation of terms due to the configuration de- 
pendence of Bf^J](RS) - once again we do not trust the 
above calculation to give the configuration dependence 
correctly. The worst contamination found for any of the 
3 samples and 3 models is \AA\ < 1.1 x 10"^. While 
the calculation presented is obviously very rough, these 
numbers do establish that the ISW/Rees-Sciama effect is 
not a significant contaminant. 



APPENDIX D: CORRELATIONS BETWEEN 
DIFFERENT FOREGROUNDS 

In the ISW analysis, for which we measure Cf"^, the 
foregrounds added linearly in the sense that Cf^ is the 
sum of the "signal" (ISW) plus the contribution from 
point sources, plus the contribution from Galactic emis- 
sion, etc. Because lensing depends on the bispectrum 
Bfvi^' foregrounds no longer add linearly: in addi- 
tion to lensing, LSS-point source-point source, and LSS- 
Galactic-Galactic terms, it is possible to have cross-terms 
involving multiple foregrounds. Out of the 6 foregrounds 
we have considered - Galactic emission, radio sources, 
infrared sources, tSZ, kSZ, and RS - it is possible to 
construct 6 x 5/2 = 15 cross-bispectra of the type LSS- 
foreground #l-foreground #2. Not all 10 combinations 
are possible: physically the foreground structure of the 
Galaxy cannot correlate with extragalactic foregrounds, 
and kSZ exhibits a reversal of sign depending on the ra- 
dial velocity that is not exhibited by point sources, tSZ, 
or RS, so correlations involving these two foregrounds 
are not possible. This leaves us with 6 possible fore- 
ground correlations involving radio point sources (rps), 
infrared point sources (irps), tSZ, and RS. As we have 
seen the RS effect is small and we will not consider its 
cross-correlations with other foregrounds here. 



1. Galaxy-radio source-infrared source 

The correlated-foreground signals involving point 
sources are most easily investigated in the context of 
the halo model. We investigate it first at a single fre- 
quency and then consider what happens in a multifre- 
quency analysis such as ours. 

The 1- and 2a-halo terms for the total foreground con- 
tain a factor of in Eq. ([55)1 where F is the flux from 
that source. In the presence of multiple emission compo- 
nents (rps, irps, tSZ) these terms contain the square of 
the total flux, i.e. 



F-^ 



(Dl) 



Therefore the contamination from galaxy-rps-irps is ob- 
tained from Eq. ([M)) by the replacement F^ 2Fi.psi^irps. 
Now the l+2a-halo corrections are all semipositive- 
definite linear combinations of F'^'s, i.e. we may schemat- 
ically write 



AA 



(D2) 



where > and the sum is over halo types fi. (The 
— sign is introduced since the contributions to are 
negative, as we are using only multipoles with the nega- 
tive sign of rps{l)/Ri and assuming positively correlated 
haloes. Also in the context of the usual halo model the 
sum would actually be an integral over the halo mass, as 
well as any additional relevant parameters; this does not 
affect the validity of the following argument.) However 
a trivial application of Cauchy's inequality (see Eq. 3.2.9 
of Ref . JTlJ) implies that 



< 2 



\ 



(D3) 

The qualitative conclusion is that the cross-talk between 
radio and infrared sources must be negligible if both types 
of sources are individually negligible; the quantitative re- 
sult is that - since we only use multipoles with the same 
sign of rps{l)/Ri ~ Eqs. ^ and jDSl imply 



I Av4 1 rps X irps (1,2a) < 2^|AA|rps(l,2a)|AA|irps(l,2a). 

(D4) 

(We have written the l-|-2a halo term with 1,2a in paren- 
theses instead of H-2ah to save space.) 

The derivation of Eq. (|D4p is valid only at one fre- 
quency, i.e. for contamination to the galaxy- VV or 
galaxy- WW correlations. If one makes no assumptions 
about the SED of the sources, then it does not need 
to be valid for the galaxy- VW correlation. [A simple 
counterexample is that if one decided to make the V- 
band flux of all infrared sources go to zero while keeping 
their W-band luminosity fixed, then the left-hand side of 
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Eq. (|D4p stays finite because of the galaxy-radio source 
(V band)-infrared source (W band) term, but the right 
side becomes zero since there is no infrared source con- 
tamination to the galaxy-Ty-Tw bispectrum.] However 
because our weights for the "TT" frequency combina- 
tion avv ■ avw ■ aww are approximately in the ratio 
0.532 . 2(0.53)(0.47) : 0.472, we should be able to approx- 
imate the galaxy-TT contamination by the replacing F 
in the above equations with 0.53Fy -I- QAlFw- The con- 
straints such as Eq. (|D4[) are then valid for the TT case. 



2. Galaxy-point source-tSZ 

Essentially the same arguments of the previous section 
can be made to the galaxy-point source-tSZ correlation. 
The one exception is that the tSZ flux is negative whereas 
the point source flux is positive. Thus for the l+2a-halo 
term Eq. (jD4p becomes 



A^|pVxtsz(l,2a) < 2J|AA|^/(l,2a)|AA|rs^(l,2 



< |A<3^(l,2a)-H|AA|rs^(l,2a). 



(D5) 



In this case, because AApsxtsz bas the opposite sign as 
AAps or Atsz (it is positive instead of negative) it fol- 
lows from this that the galaxy-ps-tSZ contamination has 
the opposite sign as the galaxy-ps-ps and galaxy-tSZ- 
tSZ contamination, but it must be smaller in magnitude 
than their sum. In other words, inclusion of the l-|-2a- 
halo galaxy-ps-tSZ term can only reduce the foreground 
contamination, and we can be conservative by ignoring 
it. 

Once again, a similar argument applies to the 2b- and 
3-halo terms, and the result should be valid for TT as 
well as for VV or WW. 



APPENDIX E: LIKELIHOOD FUNCTION 

The weak lensing of the CMB does not yet provide 
a competitive cosmological constraint. Nevertheless we 
have included it in some of the Markov chains in Paper 
I, both for completeness and as a proof of principle. This 
section describes the lensing likelihood function. 

We have included only the amplitude information from 
the galaxy-convergence cross-spectrum Cf^ (i.e. not the 
shape) since the signal-to-noise is not high enough to re- 
liably determine the latter. The amplitude is given by 
three numbers A(LRG), A(QSO), and A(NVSS) mea- 
sured in this paper, each of which is the amplitude for 
one of the samples normalized to A = \ for the fiducial 
WMAP cosmology. In order to construct a likelihood 
function, we need (i) the best values of each where 
H e{LRG,QSO,NVSS}; (ii) their covariance matrix C^"'] 
and (iii) a function to compute the expected value of 
for any choice of cosmological parameters p. Then we 



have a function 

= [C-\.K(obs) - (A^)(p)][A''(obs) - (A^)(p)]. 

(El) 

Items (i) and (ii) are easy to come by; the observed val- 
ues of A are 0.727, 1.200, and 1.110, and the inverse- 
covariance matrix is 



1.951 -0.263 -0.900 
C"^ = I -0.263 1.871 0.110 
-0.900 0.110 4.025 



> . 



(E2) 



Item (iii) is slightly harder. For each cosmology and 
each sample we can calculate the cross-power spectrum 
Cf" just as was done for the fiducial model in Eq. ([3]). 
After this we must determine the expectation values of 
the cross-spectrum estimators = Cf^ constructed in 
Section flVB II We can construct this by recalling that 
c'^ = [F^^]^^ qb where F is the response matrix and q 
is the vector of quadratic estimators. Then 



I 

I 



(E3) 



where g is the length- iVpix,LSS vector of galaxy overden- 
sities in each pixel and v is the length-2A'^pix,cMB vector 
corresponding to the reconstructed lensing map (v is a 
vector field on the sphere so there are 2 components per 
pixel). Then we have 



{QA) 



^ i?iCf "Tr [w(9)PAw(")nf 



(E4) 



It follows that the expectation values of the cross- 
spectrum estimators are 



where the quadratic estimator window function is 



w(9)p«w(^)n;^ 



From Eq. we then see that 
(^>(P 



(E5) 



(E6) 



(E7) 



where C;^^(fid) is the theoretical cross-spectrum for the 
fiducial cosmology. In this equation, only Cf,'^ needs to 
be recomputed for each new cosmology; the remaining 
coefficients can be computed once and saved. 
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